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Abstract 


A  modified  k-e  model  is  proposed  for  the  simulation  of  the  mean  wind  speed  and  turbulence 
for  a  neutrally-stratified  flow  through  and  over  a  building  array,  where  groups  of  buildings 
in  the  array  are  aggregated  and  treated  as  a  porous  barrier.  This  model  is  based  on  time 
averaging  the  spatially-averaged  Navier-Stokes  equation,  in  which  the  effects  of  the  obstacle- 
atmosphere  interaction  are  included  through  the  introduction  of  a  volumetric  momentum 
sink  (representing  drag  on  the  unresolved  buildings  in  the  array).  In  addition,  closure  of 
the  time- averaged,  spatially-averaged  Navier-Stokes  equations  requires  two  additional  prog¬ 
nostic  equations,  namely  one  for  the  time-averaged  sub-filter  kinetic  energy,  7c,  and  another 
for  the  dissipation  rate,  e,  of  7c.  The  transport  equation  for  7c  can  be  derived  from  first 
principles  and  explicitly  includes  additional  sources  and  sinks  that  arise  from  time  aver¬ 
aging  the  product  of  the  spatially-averaged  velocity  fluctuations  and  the  distributed  drag 
force  fluctuations.  The  latter  time-averaged  product  can  be  approximated  systematically 
to  any  degree  of  accuracy  using  a  Taylor  series  expansion  and,  to  this  end,  a  high-order 
approximation  is  derived  to  represent  this  source/sink  term  in  the  transport  equation  for  7c 
which  corresponds  physically  to  the  work  done  against  pressure  (form)  and  viscous  drag  in 
the  building  array.  The  dissipation  rate  U-)  equation  is  simply  obtained  as  a  dimensionally 
consistent  analog  of  the  7c-equation. 

Because  measurements  of  the  spatially-averaged  velocity  statistics  in  obstacle  arrays  are 
not  available,  the  performance  of  the  proposed  model  and  some  simplified  versions  derived 
from  it  are  compared  with  the  spatially-averaged,  time-mean  velocity  and  various  spatially- 
averaged  Reynolds  stresses  diagnosed  from  high-resolution  computational  fluid  dynamics 
(CFD)  simulations  of  the  flow  within  and  over  an  aligned  array  of  sharp-edged  cubes  with 
a  plan  area  density  of  0.25.  However,  it  should  be  emphasized  that  the  high-resolution 
CFD  flow  simulations  have  been  validated  with  wind  tunnel  experiments,  and  after  these 
validations  the  model  can  be  used  to  diagnose  the  spatially- averaged  velocity  statistics 
required  for  the  validation  of  the  distributed  drag  force  model.  It  was  found  that  the 
model  predictions  for  mean  wind  speed  and  turbulence  in  the  building  array  were  not 
sensitive  to  the  differing  treatments  of  the  source  and  sink  terms  in  the  7c-  and  e-equations, 
implying  that  the  high-order  approximations  of  these  source/sink  terms  did  not  offer  any 
predictive  advantage.  A  possible  explanation  for  this  is  the  utilization  of  the  Boussinesq 
linear  stress-strain  constitutive  relation  within  the  k-e  modelling  framework,  whose  implicit 
omission  of  any  anisotropic  eddy-viscosity  effects  renders  it  incapable  of  predicting  any 
strong  anisotropy  of  the  turbulence  field  that  might  exist  in  the  building  array.  Four  different 
methods  for  diagnosis  of  the  drag  coefficient  Co  for  the  aligned  cube  array,  required  for  the 
volumetric  drag  force  representation  of  the  cubes,  are  investigated  here. 
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Resume 


Un  modele  modifie  k-e  a  ete  propose  pour  simuler  la  vitesse  et  la  turbulence  moyennes  du 
vent  pour  des  ecoulements  d’air  de  stratification  neutre,  a  travers  et  par-dessus  une  matrice  de 
batiments.  A  I’interieur  de  la  matrice,  des  groupes  de  batiments  ont  ete  agreges  et  traites 
comme  une  barriere  poreuse.  Ce  modele  est  base  sur  le  calcul  des  moyennes  en  temps  et  en 
espace  de  liquation  Navier-Stokes  dans  lequel  les  effets  de  Pinteraction  entre  les  obstacles  et 
l’atmosph£re  sont  inclus  au  moyen  de  l’introduction  d’un  accumulateur  des  impulsions 
volumetriquese  (representant  la  trainee  sur  les  batiments  non  resolus  dans  la  matrice).  De 
plus,  la  fermeture  des  Equations  Navier-Stokes,  dont  les  moyennes  en  temps  et  espace  ont  £te 
calculees  exigent  deux  Equations  de  provision  supplementaires,  a  savoir,  une  Equation  pour 
Penergie  cinetique  sous-filtre  K  dont  les  moyennes  ont  ete  calculees  et  une  autre  Equation 
pour  le  taux  de  dissipation  s,  de  K  .  Liquation  de  transport  pour  K  peut  etre  derivee  a  partir 
des  principes  de  base  et  peut  inclure  explicitement  des  sources  et  des  puits  additionnels 
provenant  du  calcul  des  moyennes  en  temps  du  produit  des  fluctuations  des  vitesses  moyennes 
en  espace  et  des  fluctuations  de  la  force  distribute  de  trainee.  Ce  dernier  produit  aux 
moyennes  calculees  peut  etre  tvalue  systematiquement  a  n’importe  quel  niveau  d  exactitude 
en  utilisant  une  expansion  de  la  serie  de  Taylor  et,  a  cet  effet,  une  approximation  d  ordre 
superieur  est  derivee  pour  representer  ce  terme  de  la  source  ou  de  puits  dans  1  equation  de 
transport  pour  K  qui  correspond  au  travail  physiquement  accompli  contre  la  pression  (forme) 
et  la  trainee  visqueuse,  a  l’interieur  de  la  matrice  de  batiments.  L  equation  du  taux  de 
dissipation  (e-)  est  simplement  obtenue  comme  un  analogue  de  dimension  constante  de 
1 ’equation  K  . 

La  prise  de  mesures  des  statistiques  de  vitesse  moyenne  en  espace  n’etant  pas  possible  a 
l’interieur  des  matrices  d’obstacles,  le  rendement  du  modele  propose  et  de  celui  de  quelques 
versions  simplifies  qui  en  sont  derivees,  sont  compares  avec  la  moyenne  des  vitesses 
moyennes  en  espace  et  de  certains  efforts  de  Reynold  moyens  en  espace  diagnostiques  a  partir 
de  simulations  h  haute  resolution  de  la  dynamique  numerique  des  fluides  de  1  ecoulement  a 
l’interieur  et  par-dessus  la  matrice  de  cubes  a  angles  vifs  alignes  d’une  masse  surfacique  de 
0,25.  II  faut  souligner,  cependant,  que  les  simulations  d’ecoulement  de  la  dynamique 
numerique  des  fluides  h  haute  resolution  ont  ete  validees  par  des  essais  dans  des  tunnels 
a£rodynamiques  et  qu’apres  avoir  ete  valide,  le  module  peut  etre  utilise  pour  diagnostiquer  les 
statistiques  de  vitesse  moyennes  en  espacerequises  pour  la  validation  du  module  de  la  force 
distribute  de  la  trainee.  On  a  trouve  que  les  predictions  de  modeles  pour  la  vitesse  et  la 
turbulence  moyennes  du  vent  dans  la  matrice  des  batiments  ne  sont  pas^sensibles  aux 
differents  traitements  des  termes  de  la  source  et  de  puits  dans  les  equations  K  -  et  e-,  ce  qui 
implique  que  les  approximations  superieures  des  termes  de  sources  de  puits  n  offfent  aucun 
avantage  en  ce  qui  concerne  la  prediction.  II  est  possible  d’expliquer  ceci  par  la  relation 
constitutive  de  la  theorie  de  Boussinesq,  a  l’interieur  du  contexte  de  modelisation  k-t,  dont 
l’omission  implicite  des  effets  anisotropiques  de  viscosite  turbulente  le  rende  incapable  de 
predire  aucune  anisotropie  forte  du  champ  de  la  turbulence  qui  pourrait  exister  dans  la  matrice 
des  batiments.  On  etudie  ici  quatre  methodes  differentes  de  diagnostic  du  coefficient  Qj  de  la 
trainee  concemant  la  matrice  de  cubes  alignes  requise  pour  representer  la  force  de  la  trainee 
des  cubes. 
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Executive  summary 


Introduction:  It  is  anticipated  that  Canadian  Forces  (CF)  in  the  foreseeable  future  will 
have  to  fight  in  or  protect  urban  areas,  whether  in  battle,  peace-making,  peacekeeping, 
or  counter-terrorist  operations.  The  increased  awareness  and  importance  accorded  by  the 
public  worldwide  and  their  governments  to  maintain  appropriate  defences  against  chemical 
and  biological  warfare  (CBW)  agents  in  an  urban  (built-up)  environment,  the  prediction 
of  casualties  and  human  performance  degradation  resulting  from  such  releases,  and  the 
development  of  operational  procedures  and  regulations  to  control,  mitigate,  and  monitor  the 
fate  of  CBW  agents  in  urban  areas  with  high  population  densities,  will  require  mathematical 
modeling  of  urban  wind  flows  and  dispersion.  In  this  regard,  it  should  be  noted  that  the 
prediction  of  flows  in  the  urban  environment  is  in  principle  pre-requisite  to  or  co-requisite 
with  the  prediction  of  contaminant  (e.g.,  CBW  agent)  dispersion  within  a  cityscape. 

Results:  A  modified  k-e  model  is  proposed  for  the  simulation  of  the  mean  wind  speed 
and  turbulence  for  a  neutrally-stratified  flow  through  and  over  a  building  array,  where 
groups  of  buildings  in  the  array  are  aggregated  and  treated  as  a  porous  barrier.  This 
model  is  based  on  time  averaging  the  spatially-averaged  Navier-Stokes  equation,  in  which 
the  effects  of  the  obstacle-atmosphere  interaction  are  included  through  the  introduction  of  a 
volumetric  momentum  sink  (representing  drag  on  the  unresolved  buildings  in  the  array).  In 
addition,  closure  of  the  time-averaged,  spatially-averaged  Navier-Stokes  equations  requires 
two  additional  prognostic  equations,  namely  one  for  the  time-averaged  sub-filter  kinetic 
energy,  k,  and  another  for  the  dissipation  rate,  e,  of  k.  The  transport  equation  for  7c  can  be 
derived  from  first  principles  and  explicitly  includes  additional  sources  and  sinks  that  arise 
from  time  averaging  the  product  of  the  spatially-averaged  velocity  fluctuations  and  the 
distributed  drag  force  fluctuations.  The  latter  time-averaged  product  can  be  approximated 
systematically  to  any  degree  of  accuracy  using  a  Taylor  series  expansion  and,  to  this  end, 
a  high-order  approximation  is  derived  to  represent  this  source/sink  term  in  the  transport 
equation  for  k  which  corresponds  physically  to  the  work  done  against  pressure  (form)  and 
viscous  drag  in  the  building  array.  The  dissipation  rate  (e-)  equation  is  simply  obtained  as 
a  dimensionally  consistent  analog  of  the  ^-equation. 

The  performance  of  the  proposed  model  and  some  simplified  versions  derived  from  it  are 
compared  with  the  spatially-averaged,  time-mean  velocity  and  various  spatially-averaged 
Reynolds  stresses  diagnosed  from  high- resolution  computational  fluid  dynamics  (CFD)  sim¬ 
ulations  of  the  flow  within  and  over  an  aligned  array  of  sharp-edged  cubes  with  a  plan  area 
density  of  0.25.  It  was  found  that  the  model  predictions  for  mean  wind  speed  and  turbu¬ 
lence  in  the  building  array  were  not  sensitive  to  the  differing  treatments  of  the  source  and 
sink  terms  in  the  7c-  and  e-equations,  implying  that  the  high-order  approximations  of  these 
source/sink  terms  did  not  offer  any  predictive  advantage.  A  possible  explanation  for  this 
is  the  utilization  of  the  Boussinesq  linear  stress-strain  constitutive  relation  within  the  k-e 
modelling  framework,  whose  implicit  omission  of  any  anisotropic  eddy-viscosity  effects  ren¬ 
ders  it  incapable  of  predicting  any  strong  anisotropy  of  the  turbulence  field  that  might  exist 
in  the  building  array.  Four  different  methods  for  diagnosis  of  the  drag  coefficient  Co  for 
the  aligned  cube  array,  required  for  the  volumetric  drag  force  representation  of  the  cubes, 
are  investigated  here. 
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Significance  and  Future  Plans:  A  knowledge  of  the  structure  of  the  mean  flow  and 
turbulence  describing  the  complex  flow  patterns  within  and  over  clusters  of  buildings  is 
essential  for  improving  urban  dispersion  models.  Unfortunately,  the  computational  demands 
of  computational  fluid  dynamics  (CFD)  where  all  buildings  are  resolved  explicitly  in  the 
sense  that  boundary  conditions  are  imposed  at  all  surfaces  (e.g.,  walls,  roofs)  are  so 
prohibitive  as  to  preclude  their  use  for  emergency  response  situations  which  require  the  ability 
to  generate  an  urban  flow  and  dispersion  prediction  in  a  time  frame  that  will  permit  protective 
actions  to  be  taken.  In  view  of  this,  we  have  demonstrated  in  this  report  how  to  construct 
models  for  the  statistics  of  the  mean  flow  and  turbulence  in  an  urban  canopy  that  are  obtained 
by  averaging  horizontally  the  mean  wind  and  turbulence  statistics  over  an  area  that  is  larger 
than  the  spacings  between  the  individual  roughness  elements  comprising  the  urban  canopy, 
but  less  than  the  length  scales  over  which  the  roughness  element  density  changes.  The 
development  of  spatially-averaged  Reynolds-averaged  Navier-Stokes  models  permits  an 
efficient  prediction  of  urban  flows  required  for  emergency  response  situations.  The  utility  of 
the  simplified  urban  flow  simulation  models  investigated  here  for  provision  of  the  disturbed 
wind  field  statistics  required  by  a  physically-based  dispersion  model  needs  further 
investigation. 


Yee  E  and  Lien,  F.S.  (2004).  A  Distributed  Drag  Force  Approach  for  the  Numerical 
Simulation  of  Urban  Flows.  (DRDC  Suffield  TR  2004-169).  Defence  R&D  Canada  -  Suffield. 
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Sommaire 


Introduction  :  On  prevoit  que  les  Forces  canadiennes  (FC),  auront,  dans  un  avenir  assez 
rapproche,  a  combattre  a  l’interieur  de  zones  urbaines  ou  a  proteger  ces  dernieres,  durant  des 
operations  de  combat,  de  maintien  de  la  paix  ou  antiterroristes.  Une  prise  de  conscience 
accrue  et  l’importance  que  le  monde  entier  et  ses  gouvernements  accordent  au  maintien  de 
moyens  de  defense  appropries  contre  les  agents  de  guerre  chimiques  et  biologiques  (CB)  dans 
les  milieux  urbains  (construits),  a  la  prevision  des  blesses  et  la  degradation  de  la  performance 
humaine  resultant  de  telles  emissions  ainsi  qu’a  la  mise  au  point  de  procedures 
op^rationnelles  et  de  reglements  de  controle  pour  attenuer  et  surveiller  le  sort  des  agents  CB 
dans  les  zones  urbaines  comprenant  des  hautes  densites  de  population,  exigeront  des 
mod^lisations  math^matiques  des  ecoulements  eoliens  urbains  et  de  leur  dispersion.  A  cet 
£gard,  il  faut  noter  que  la  prevision  des  ecoulements  dans  un  milieu  urbain  est  en  principe  pr£- 
requise  ou  co-requise  avec  la  prevision  de  la  dispersion  du  contaminant  (par  ex.  :  agent  CB), 
dans  un  paysage  urbain. 

Resultats  :  Un  modele  modifie  k- s  a  ete  propose  pour  simuler  la  vitesse  et  la  turbulence 
moyennes  du  vent  pour  des  ecoulements  d’air  de  stratification  neutre,  a  travers  et  par-dessus 
une  matrice  de  batiments.  A  l’interieur  de  la  matrice,  des  groupes  de  batiments  ont  ete  agreges 
et  traites  comme  une  barriere  poreuse.  Ce  modele  est  base  sur  le  calcul  des  moyennes  en 
temps  et  en  espace  de  l’equation  Navier-Stokes,  dans  lequel  les  effets  de  1’interaction  entre  les 
obstacles  et  l’atmosphere  sont  inclus  au  moyen  de  1’ introduction  d’un  accumulateur  des 
impulsions  volumetriques  (representant  la  trainee  sur  les  batiments  non  resolus  dans  la 
matrice).  De  plus,  la  fermeture  des  equations  Navier-Stokes,  dont  les  moyennes  en  temps  et 
en  espace  ont  et6  calculees,  exigent  deux  Equations  de  prevision  supplementaires,  a  savoir, 
une  Equation  pour  l’energie  cinetique  sous-filtre  K  dont  les  moyennes  ont  ete  calculees  et  une 
autre  equation  pour  le  taux  de  dissipation  e,  de  K  .  L’equation  de  transport  pour  K  peut  etre 
derivee  a  partir  des  premiers  principes  et  peut  inclure  explicitement  des  sources  et  des  puits 
additionnels  provenant  du  calcul  des  moyennes  en  temps  du  produit  des  fluctuations  des 
vitesses  moyennes  en  espace  et  des  fluctuations  de  la  force  distribute  de  trainee.  Ce  dernier 
produit  aux  moyennes  calculees  peut  etre  evalue  systematiquement  a  n’importe  quel  niveau 
d’exactitude  en  utilisant  une  expansion  de  la  serie  de  Taylor  et,  a  cet  effet,  une  approximation 
d’ordre  suptrieur  est  derivee  pour  representer  ce  terme  de  la  source  ou  de  piegeage  dans 
1’equation  de  transport  pour  K  qui  correspond  au  travail  physiquement  accompli  contre  la 
pression  (forme)  et  la  trainee  visqueuse,  a  l’interieur  de  la  matrice  de  batiments.  L’equation  du 
taux  de  dissipation  (e-)  est  obtenue  simplement  comme  un  analogue  de  dimension  constante 
de  l’equation  K . 

Le  rendement  du  modele  propose  et  de  quelques  versions  simplifies  qui  en  sont  d^rivees, 
sont  compares  entre  la  moyenne  des  vitesses  instantanees,  dont  les  moyennes  en  espace  ont 
et6  calculees  et  des  tensions  de  Reynold  variees  dont  les  moyennes  en  espace  ont  ete 
calculees,  diagnostiquees  a  partir  de  simulations  a  haute  resolution  de  la  dynamique 
numerique  des  fluides  de  I’ecoulement  a  l’interieur  et  par-dessus  la  matrice  de  cubes  a  angles 
vifs  alignes  d’une  masse  surfacique  de  0,25.  On  a  trouve  que  les  predictions  de  modeles  pour 
la  vitesse  et  la  turbulence  moyennes  du  vent  dans  la  matrice  des  batiments  ne  sont  pas 
sensibles  aux  differents  traitements  des  termes  de  la  source  et  de  puits  dans  les  equations  K  - 
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et  £-,  ce  qui  implique  que  les  approximations  superieures  des  termes  de  sources  /de  puits 
n’offfent  aucun  avantage  en  ce  qui  conceme  la  prediction.  II  est  possible  d’expliquer  ceci  par 
la  relation  constitutive  de  la  theorie  de  Boussinesq,  a  Pinterieur  du  contexte  de  modelisation 
k-E,  dont  l’omission  implicite  des  effets  anisotropiques  de  viscosite  turbulente  le  rende 
incapable  de  predire  aucune  anisotropie  forte  du  champ  de  la  turbulence  qui  pourrait  exister 
dans  la  matrice  des  batiments.  On  etudie  ici  quatre  methodes  differentes  de  diagnostic  du 
coefficient  Cq  de  la  trainee  pour  la  matrice  de  cubes  alignes  qui  est  requise  pour  reprdsenter  la 

force  de  la  trainee  des  cubes. 

La  portee  des  resultats  et  les  plans  futurs  :  II  est  essentiel  de  connaltre  la  structure  de 
Pecoulement  moyen  et  de  la  turbulence  moyenne  qui  decrit  les  modeles  complexes 
d’^coulement  a  Pinterieur  et  par-dessus  des  r^seaux  de  batiments  et  qui  vise  a  amdliorer  les 
modules  urbains  de  dispersion.  Malheureusement,  tous  les  batiments  etant  resolus 
explicitement  avec  des  conditions  de  limites  imposees  sur  toutes  les  surfaces  (par  ex. .  murs, 
toits),  les  exigences  en  calculs  de  la  dynamique  numerique  des  fluides  rendent  impossible 
Putilisation  de  cette  methode  en  cas  d’une  intervention  d’urgence  ;  une  situation 
d’intervention  d’urgence  exige  de  generer  une  prediction  des  ecoulements  et  de  la  dispersion 
en  milieu  urbain  dans  un  delai  d’execution  assez  bref  pour  permettre  la  prise  des  mesures  de 
protection.  Dans  cette  optique,  nous  avons  d^montre  dans  ce  rapport,  comment  construire  des 
modeles  de  statistiques  concemant  Pecoulement  moyen  et  les  turbulences  moyennes  dans  une 
couverture  urbaine.  Ces  modeles  sont  obtenus  en  calculant  la  moyenne  horizontale  des 
statistiques  de  Pecoulement  moyen  et  des  turbulences  moyennes  du  vent,  sur  un  territoire  plus 
large  que  Pespacement  entre  les  elements  individuels  de  rugosite  comprenant  la  couverture 
urbaine  mais  moins  large  que  les  echelles  de  longueur  couvrant  un  territoire  sur  lequel  change 
la  densite  de  Peidment  de  rugosite.  La  mise  au  point  des  modeles  Reynold  et  Navier-Stokes, 
dont  les  moyennes  en  espace  ont  ete  calcuiees,  permet  une  prediction  efficace  des 
ecoulements  urbains  pour  des  situations  d’intervention  d’urgence.  L’utilite  des  modeles 
simplifies  de  simulation  de  Pecoulement  urbain  ayant  ete  etudiee  ici  pour  les  statistiques  de 
champs  de  vent  perturbe  requis  par  le  modele  de  dispersion,  d’apres  des  criteres  physiques, 
exige  d’etre  etudiee  plus  profondement. 


Yee,  E.  and  Lien,  F.S.  (2004).  A  Distributed  Drag  Force  Approach  for  the  Numerical  Simulation 
of  Urban  Flows.  (DRDC  Suffield  TR  2004-169).  R  &  D  pour  la  defense  Canada  -  Suffield. 
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Introduction 


The  turbulent  flow  within  and  over  urban  areas  covered  with  agglomerations  of  discrete 
buildings,  often  with  irregular  geometry  and  spacing,  is  generally  very  complex  and  pos¬ 
sesses  a  fully  three-dimensional  structure.  Although  the  application  of  computational  fluid 
dynamics  (CFD)  to  the  prediction  of  the  mean  flow  and  turbulence  near  and  around  a  single 
building  or  within  and  over  a  regular  array  (or,  canopy)  of  buildings  is  progressing  ([1],  [2], 
[3],  [4]),  this  method  tends  to  require  extensive  computational  resources.  Nevertheless,  CFD 
simulations  which  involve  the  solution  of  the  conservation  equations  for  mass,  momentum, 
and  energy  allow  the  prognosis  of  a  number  of  velocity  statistics  (e.g.,  mean  velocity,  nor¬ 
mal  stresses,  shear  stresses,  etc.)  in  an  urban  canopy.  A  knowledge  of  the  structure  of  the 
mean  flow  and  turbulence  describing  the  complex  flow  patterns  within  and  over  clusters  of 
buildings  is  also  essential  for  improving  urban  dispersion  models. 

Unfortunately,  the  computational  demands  of  CFD  where  all  buildings  are  resolved  expli¬ 
citly  in  the  sense  that  boundary  conditions  are  imposed  at  all  surfaces  (e.g.,  walls,  roofs) 
are  so  prohibitive  as  to  preclude  their  use  for  emergency  response  situations  which  require 
the  ability  to  generate  an  urban  flow  and  dispersion  prediction  in  a  time  frame  that  will 
permit  protective  actions  to  be  taken.  In  view  of  this,  we  argue  that  for  many  practical 
applications  it  is  convenient  to  consider  the  prediction  of  spatial  averages  of  the  mean  wind 
and  turbulence  in  an  urban  canopy  that  is  obtained  by  averaging  horizontally  the  mean 
wind  and  turbulence  statistics  over  an  area  that  is  larger  than  the  spacings  between  the 
individual  roughness  elements  comprising  the  urban  canopy,  but  less  than  the  length  scale 
over  which  the  roughness  element  density  changes. 

In  this  report,  we  focus  on  the  mathematical  formulation  of  a  numerical  model  for  the 
prediction  of  flows  within  and  over  a  building  array  based  on  an  aggregation  of  groups  of 
buildings  in  the  array  into  a  number  of  ‘drag  units’,  with  the  ensemble  of  units  being  treated 
as  a  continuous  porous  medium.  This  approach  will  obviate  the  need  to  impose  boundary 
conditions  along  the  surfaces  of  all  buildings  (and  other  obstacles)  in  the  array.  Wilson 
and  Yee  [5]  applied  something  like  this  approach  to  simulate  the  mean  wind  and  turbulence 
energy  fields  in  a  single  unit  cell  of  the  wind  tunnel  “Tombstone  Canopy”  [6],  a  regular 
diamond  staggered  array  of  bluff  (impermeable)  aluminum  plates,  with  a  disappointing 
outcome  (subsequent  work  showed  that  invoking  a  Reynolds  stress  closure  did  not  help) .  We 
now  know  this  may  owe  to  the  existence  of  (previously  unsuspected)  large  eddies  generated 
by  the  strong  shear  layer  near  the  top  of  the  canopy,  eddies  that  span  more  than  one  unit 
cell  in  the  streamwise  direction,  and  imply  that  imposition  of  an  artificial  condition  of 
periodicity  at  the  boundaries  of  a  single  cell  amounts  to  solving  a  different  flow  problem. 
These  large-scale  (coherent)  eddy  structures  generated  at  or  near  the  canopy  top  have  been 
observed  using  highly  resolved,  two-dimensional  laser-induced  fluorescence  measurements  of 
the  fine  structure  of  the  fully  space-  and  time- varying  conserved  scalar  field  resulting  from  a 
point-source  release  of  a  tracer  within  the  “Tombstone  Reloaded  Canopy”  in  a  water  channel 
simulation  [7].  Belcher  et  al.  [8]  applied  a  similar  approach  to  investigate  the  adjustment  of 
the  mean  velocity  to  a  canopy  of  roughness  elements  using  a  linearized  flow  model  (obtained 
by  determining  analytically  small  perturbations  to  the  undisturbed  upstream  logarithmic 
mean  velocity  profile  induced  by  the  drag  due  to  an  obstacle  array).  Hookham  et  al.  [9] 
have  used  a  drag  force  parameterization  to  represent  the  effects  of  buildings  on  the  flow  in 
the  development  of  their  Urban  Windfield  Module. 
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There  is  precedent  for  treating  drag  on  unresolved  buildings  in  an  urban  canopy  by  means 
of  a  distributed  momentum  sink  for  the  representation  of  the  effects  on  the  mean  flow  and 
turbulence  arising  from  the  form  and  viscous  drag  on  canopy  elements.  As  motivation,  we 
recall  that  a  similar  approach  has  been  applied  over  the  past  50  years  to  the  modelling  of 
flows  in  plant  canopies  and  about  porous  windbreaks.  Although  a  sink  or  drag  term  has 
been  added  in  an  ad  hoc  fashion  to  the  free-air  mean  momentum  equation  to  model  the 
canopy  mean  wind  profile  over  a  number  of  years  ([10],  [11],  [12]),  it  was  not  until  1977 
that  Wilson  and  Shaw  [13]  showed  how  to  apply  a  rigorous  spatial-averaging  procedure  to 
obtain  the  equations  for  a  spatially-continuous  area-averaged  mean  wind  and  turbulence 
field.  In  this  seminal  work,  Wilson  and  Shaw  [13]  demonstrated  how  additional  source  and 
sink  terms  representing  the  flow  interaction  with  the  canopy  elements  emerge  naturally  by 
application  of  a  particular  spatial  averaging  procedure  to  the  Reynolds-averaged  Navier- 
Stokes  equations  that  obtain  at  every  point  in  the  canopy  airspace.  This  procedure  was 
further  developed  by  Raupach  and  Shaw  [14]  for  the  case  of  a  horizontal  plane  averaging 
operation.  In  particular,  Raupach  and  Shaw  [14]  discuss  two  different  options  for  averaging 
over  a  horizontal  plane;  namely,  horizontally  averaging  the  equations  of  motion  at  a  single 
time  instant  over  a  plane  extensive  enough  “to  eliminate  variations  due  to  canopy  structure 
and  the  largest  length  scales  of  the  turbulent  flow”  (scheme  I)  and  conventional  time  aver¬ 
aging  of  the  equations  of  motion  followed  by  horizontal  averaging  over  a  plane  large  enough 
“to  eliminate  variations  in  the  canopy  structure”  (scheme  II).  Scheme  I  has  rather  limited 
applicability  since  it  cannot  be  applied  to  horizontally  inhomogeneous  canopies. 


Finnigan  [15]  and  Raupach  et  al.  [6]  investigated  the  volume-averaging  method.  Finnigan 

[15]  considered  details  such  as  plant  motion  (e.g.,  coherently  waving  plant  canopies)  which 
gives  rise  to  a  ‘waving  production’  term  in  the  transport  equations  for  turbulence  quantities. 
We  note  that  plant  motion  is  not  a  factor  directly  pertinent  to  the  present  work  which 
focusses  on  urban  canopy  flows,  but  these  concepts  may  have  a  bearing  on  the  case  of 
moving  obstacles  (e.g.,  vehicles)  within  the  urban  canopy.  Following  ideas  of  Hanjalic  et  al. 

[16]  and  parallelling  Shaw  and  Seginer  [17],  Wilson  [18]  developed  an  empirical  two-band 
model  for  the  turbulence  kinetic  energy  (TKE)  which  represented  the  large-  and  fine-scale 
components  of  the  turbulence  and  their  dynamics  [the  multiple  time-scale  approach  has 
seen  much  subsequent  use  [19],  but  parameterizing  the  exchange  of  kinetic  energy  between 
the  spectral  bands  is  a  pre-eminent  difficulty  of  the  approach] .  Here,  the  turbulence  kinetic 
energy  was  separated  into  two  wave-bands,  corresponding  to  shear  kinetic  energy  (SKE, 
low-frequency  band)  and  wake  kinetic  energy  (WKE,  high-frequency  band),  with  separate 
equations  developed  to  represent  their  dynamics.  Wilson  [20],  Green  [21],  Wang  and  Takle 
[22],  Wang  and  Takle  [23],  Liu  et  al.  [24],  Ayotte  et  al.  [25],  Sanz  [26],  and  Wilson  and  Yee 
[27]  investigated  various  modifications  of  the  k-e  model  or  the  Reynolds  stress  transport 
model  to  account  for  interaction  of  the  air  with  canopy  elements. 


This  is  the  final  report  of  two  describing  the  numerical  modelling  of  the  developing  flow 
within  and  over  a  3-D  building  array.  In  the  first  report  (henceforth  I)  [4],  we  used  the 
Reynolds-averaged  Navier-Stokes  (RANS)  equations  in  conjunction  with  a  twoequation 
turbulence  model  (i.e.,  k-e  model)  to  predict  the  complex  three-dimensional  disturbed  flow 
within  and  over  a  3-D  building  array  under  neutral  stability  conditions.  The  simulations 
of  the  mean  flow  field  and  turbulence  kinetic  energy  were  validated  with  data  obtained 
from  a  comprehensive  wind  tunnel  experiment  conducted  by  Brown  et  al.  [28].  Here, 
it  was  demonstrated  that  the  mean  flow  and  turbulence  kinetic  energy  from  the  numerical 
and  physical  simulations  exhibited  striking  resemblences.  In  addition,  the  importance  of  the 
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kinematic  ‘dispersive  stresses’  relative  to  the  spatially-averaged  kinematic  Reynolds  stresses 
for  developing  flow  within  and  over  an  urban-like  roughness  array  has  been  quantified  using 
the  high-resolution  CFD  results  obtained  with  the  high-Reynolds-number  k-e  model. 

In  this  report,  we  focus  on  the  formulation  of  a  numerical  model  for  the  prediction  of  flows 
within  and  over  a  building  array  based  on  an  aggregation  of  groups  of  buildings  in  the 
array  into  a  number  of  ‘drag  units’,  with  each  unit  being  treated  as  a  porous  barrier.  This 
approach  will  obviate  the  need  to  impose  boundary  conditions  along  the  surfaces  of  all 
buildings  (and  other  obstacles)  in  the  array.  Here,  we  present  details  of  the  mathematical 
framework  required  to  derive  the  transport  equation  for  the  time  average  of  the  locally- 
spatially-averaged  velocity  through  a  building  array  (which  is  treated  here  as  a  porous 
medium),  and  the  two  additional  prognostic  equations  required  to  close  this  equation  set. 
These  additional  equations  predict  the  time-averaged  resolved-scale  kinetic  energy  of  turbu¬ 
lence,  k,  and  its  dissipation  rate,  e.  The  closure  problem  relating  to  the  ‘correct’  represen¬ 
tation  of  the  additional  source/sink  terms  in  the  transport  equations  for  mean  momentum, 
turbulence  energy,  and  dissipation  rate  is  investigated  in  detail.  Most  of  the  work  reported 
is  motivated  by  conceptual  and  logical  difficulties  in  the  self-consistent  treatment  of  source 
and  sink  terms  in  the  transport  equations  for  turbulence  kinetic  energy  and  its  dissipation 
rate.  To  this  end,  we  attempt  to  lay  the  foundations  for  a  systematic  mathematical  for¬ 
mulation  that  could  be  used  to  construct  the  additional  source/sink  terms  in  the  transport 
equations  for  7c  and  e,  in  response  (and  to  some  extent,  contradiction)  to  the  assertions 
made  by  Wilson  and  Mooney  [29]  that  it  is  “impossible  to  know  the  ‘correct’  influence  of 
the  unresolved  processes  at  the  fence  on  TKE  and  its  dissipation  rate”  and  by  Wilson  et  al. 
[30]  that  “ k-e  closures  give  predictions  that  are  sensitive  to  details  of  ambiguous  choices”. 
The  problem  of  the  diagnosis  of  the  drag  coefficient  Cp  required  in  the  parameterization 
of  the  distributed  drag  is  addressed.  We  show  detailed  comparisons  of  predictions  obtained 
from  the  model  with  spatial  averages  of  the  mean  velocity  field  and  second-order  velocity 
statistics  obtained  from  a  high-resolution  CFD  simulation  of  flow  within  and  over  a  3-D 
building  array. 

Model  formulation 


Spatial  and  time  averaging  operations 

Before  we  begin,  we  present  a  short  note  on  the  notation  that  will  be  used.  The  following 
derivations  will  invariably  use  the  flexibility  of  the  Cartesian  tensor  notation,  with  Roman 
indices  such  as  i,  j,  or  k  taking  values  of  1,  2,  or  3.  We  shall  also  employ  the  Einstein 
summation  convention  in  which  repeated  indices  are  summed.  For  any  flow  variable  <j b, 
{(j>)  will  denote  the  spatial  (volume)  average,  (j)  the  time  average,  <f>'  the  departure  of  (f> 
from  its  time-averaged  value,  and  <j>"  the  departure  of  <f>  from  its  spatially-averaged  value. 
In  addition,  Ui  is  the  total  velocity  in  the  Zj-direction,  with  i  —  1,  2,  or  3  representing 
the  streamwise  x,  spanwise  y,  or  vertical  z  direction.  Finally,  x  =  ( x,y,z ),  (ui,u.2,uz)  = 
( u,v,w ),  and  t  denotes  time. 

The  derivation  of  a  model  for  the  spatially-averaged  time-mean  flow  can  start  either  from 
applying  the  spatial  averaging  operation  to  the  time-averaged  Navier-Stokes  (NS)  equation 
((NS)),  or  the  time  averaging  operation  to  the  spatially-averaged  Navier-Stokes  equation 
( (NS) ).  In  formulating  such  equations,  we  must  choose  a  suitable  decomposition  of  a 
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flow  property  into  its  rapidly  and  slowly  varying  components  and  determine  a  strategy  for 
applying  the  corresponding  averaging  operation.  First,  consider  spatial  averaging  in  some 
multiply  connected  space.  In  a  “slow  +  fast”  decomposition  of  a  flow  property  4>  based  on 
spatial  filtering,  scales  are  separated  by  applying  a  low-pass  scale  filter  to  give  a  filtered 
quantity  { <fi )  defined  by 


{4>)(x)=  f  G(x  -y)4>(y)dy  =  G*4>.  (1) 

J  a.s. 

The  integral  in  Equation  (1)  is  assumed  to  be  over  all  space  (a.s.).  Here,  G(x  —  y),  the 
convolution  filter  kernel,  is  a  localized  function  (i.e.,  G  — *  0  as  ||x  —  y||  — >  oo,  where  ||  •  || 
denotes  the  Euclidean  norm)  with  width  A  which  is  related  to  some  cutoff  scale  in  space, 
and  *  is  used  to  denote  the  convolution  operation.  In  general,  the  filter  width  can  depend 
on  x,  which  we  will  explicitly  indicate  using  the  notation  G(x  —  y|A(x)). 


If  we  assume  that  G  is  a  symmetric  function  of  x  —  y,  and  differentiate  Equation  (1)  with 
respect  to  Xi,  we  get  the  following  relationship  between  the  spatial  average  of  the  spatial 
derivative  and  the  spatial  derivative  of  the  spatial  average  of  a  quantity  4>  (on  application 
of  the  Gauss  divergence  theorem): 


d(4>)  s 

dxi 


d 


dxi 

(QG  A 


<t> 


) 


da 

dxi 


+ 


J  G(x-y|A(x))0(y)ni  dS, 


(2) 


where  S  denotes  the  sum  of  all  obstacle  surfaces  contained  in  the  multiply  connected  region 
(extending  over  all  space),  zij  is  the  unit  outward  normal  in  the  z-th  direction  on  the  surface 
S  (positive  when  directed  into  the  obstacle  surface),  and  [/, g]  =  f 9  ~  df  denotes  the 
commutator  bracket  of  two  operators  /  and  g.  Equation  (2)  will  be  referred  to  as  the 
generalized  spatial  averaging  theorem.  A  special  case  of  this  theorem  (known  as  the  spatial 
averaging  theorem)  has  been  derived  by  Raupach  and  Shaw  [14]  and  Howes  and  Whitaker 

[31]. 


The  spatial  filtering  operation  does  not  commute  with  spatial  differentiation.  The  non¬ 
commutation  of  these  two  operations  results  from  two  contributions.  The  first  contribution 
is  encapsulated  in  the  first  term  on  the  right-hand-side  of  Equation  (2)  which  arises  from  the 
spatial  variation  in  filter  cutoff  length.  The  second  contribution,  summarized  in  the  second 
term  on  the  right-hand-side  of  Equation  (2),  is  due  to  the  presence  of  obstacle  surfaces  in  the 
multiply  connected  flow  domain.  Interestingly,  if  we  apply  the  spatial-averaging  operator  to 
the  continuity  equation,  the  spatial  variation  in  the  filter  width  implies  that  (ui)  is  no  longer 
solenoidal.  More  specifically,  although  the  velocity  across  the  air/solid  boundaries  vanishes 
owing  to  the  no-slip  and  impermeability  boundary  conditions  here,  the  spatial  variation  of 
the  filter  width  implies  an  extra  source/sink  term  in  the  filtered  continuity  equation  [which 
is  a  direct  consequence  of  Equation  (2)]: 


d(uj) 

dxi 


Ui  = 


(3) 
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To  ensure  that  the  spatially-averaged  velocity  field  is  solenoidal,  we  consider  a  special  con¬ 
volution  kernel  whose  filter  cutoff  length  does  not  depend  on  x.  To  this  purpose,  consider 
the  box  or  top-hat  filter  defined  as 


r(*  -  v1)  =  /  1/V*  if  \xi  ~Vi\  <  Ai/2; 
'  \  0,  otherwise. 


(4) 


Here,  V  =  AxAyAz  is  the  constant  volume  over  which  we  average  to  obtain  continuous 
variables.  With  the  constant  width  filter  kernel  of  Equation  (4),  the  spatial  average  of  a 
flow  property  <p  of  Equation  (1)  becomes  simply 


(0)(x,t)  =  ^  f^cpdV  =  ^  f  4>(x  +  r,t)dr.  (5) 

Note  that  in  Equation  (5) ,  the  averaging  volume  includes  both  fluid  and  solid  parts  (obsta¬ 
cles).  Applying  the  volume-averaging  operator  of  Equation  (5)  to  the  continuity  equation 
results  in  a  spatially-averaged  velocity  field  (u{)  that  is  solenoidal. 


We  will  use  the  spatial- averaging  operation  displayed  in  Equation  (5),  where  the  average 
is  taken  over  both  the  fluid  and  solid  phases  in  V  (averaging  volume),  and  the  normalizing 
factor  is  the  total  volume  V.  For  two-phase  systems,  two  other  definitions  for  averaging  have 
been  proposed  (e.g.,  Miguel  et  al. [32]).  In  a  two-phase  system,  the  total  averaging  volume 
V  is  made  up  of  the  volume  of  the  fluid  phase  Vf  and  the  solid  phase  Vs ,  so  V'  =  Vf  +  Vs. 
The  superficial  (external)  phase  average  of  <f>  is  defined  as 

(<P)e  =  ^jv  4>dV  (6) 

and  the  intrinsic  (internal)  phase  average  of  (j>  is  defined  as 


The  intrinsic  phase  average  is  an  average  of  a  flow  property  over  the  fluid  phase  (i.e.,  the 
averaging  volume  Vf  excludes  the  solid  phase,  with  the  normalizing  factor  being  V/).  On 
the  other  hand,  the  external  phase  average  is  a  weighted  average  of  the  fluid  property  over 
the  fluid  phase  (i.e.,  the  total  volume  V  is  used  as  the  normalizing  factor,  but  the  averaging 
excludes  the  solid  phase). 


The  spatial  (or,  volume)  average  defined  in  Equation  (5)  seems  natural  in  the  present 
context,  and  leads  to  the  simplest  forms  for  the  volume- averaged  transport  equations  on 
application  of  the  volume-averaging  operator  to  the  continuity  and  the  Reynolds  equation 
for  mean  momentum  at  a  single  point.  Although  V  is  a  constant  that  is  independent  of  the 
spatial  coordinates,  Vf  which  represents  the  volume  of  the  fluid  phase  contained  within  V 
need  not  be  (e.g.,  for  an  inhomogeneous  canopy,  Vf  and  Vs  will  be  a  function  of  the  spatial 
coordinates).  Because  Vf  depends  on  x  in  Equation  (7)  (i.e.,  the  filter  width  depends 
on  x),  the  use  of  the  intrinsic  phase  average  will  result  in  a  filtered  velocity  that  is  not 
solenoidal  [cf.  Equation  (3)].  Even  so,  transport  equations  for  the  intrinsic  phase-averaged 
velocity  ( Uj)i  can  be  derived,  provided  that  the  fluid-phase  volume  V/(x)  is  differentiable 
(or,  equivalently,  that  the  porosity  £  =  Vf/V  is  a  differentiable  function  of  x)  although  this 
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will  result  in  a  number  of  ‘extra’  source/sink  terms  in  these  equations  arising  solely  from 
the  dependence  of  Vf  on  the  spatial  coordinates.  The  external  phase  average  does  not  seem 
natural  in  the  present  context  because  ((<f>)e)e  7^  (4>)el  (ori  equivalently,  (<pe)e  r1  0  where 

<t>"  =  <i>-  {4>)e)- 

In  analogy  with  the  spatial  (or,  volume)  average  defined  in  Equation  (5),  the  time  average 
of  0,  which  we  denote  using  an  overbar,  will  be  defined  as 

1  rto+T 

<£(x)  =  ^  /  4>(x,  t )  dt,  Ti  <  T «:  T2.  (8) 

Jto 

In  view  of  Equations  (5)  and  (8),  time  and  spatial  averaging  commute  so  (<j>)  =  {$).  In 
Equation  (5),  the  horizontal  averaging  scales  Ax,  Ay  need  to  be  large  compared  to  the 
separation  between  individual  roughness  elements,  but  much  less  than  the  characteristic 
length  scales  over  which  the  density  of  the  roughness  elements  changes;  but  to  ensure  a 
sufficient  vertical  resolution  of  the  flow  property  gradients,  Az  -C  Ax,Ay  making  V  a 
thin,  horizontal  slab.  In  Equation  (8),  the  averaging  time  T  is  implicitly  assumed  to  be 
sufficiently  long  to  ensure  that  many  cycles  of  the  rapid  turbulent  fluctuations  in  a  flow 
property  are  captured,  but  sufficiently  short  so  that  the  external  large-scale  variations  in 
the  flow  property  are  approximately  constant.  Hence,  in  Equation  (8),  T\  and  T2  denote 
the  time  scales  characteristic  of  the  rapid  and  slow  variations  in  the  flow  property  <f),  with 
the  implicit  assumption  that  T\  and  T2  differ  by  several  orders  of  magnitude. 

In  general,  <fi  can  be  decomposed  in  the  following  two  ways: 


(j>  =  4>  +  <t>'i  4>'  —  0> 

(9) 

*=(*)  +  *",  (</>")  =0. 

(10) 

Although  (NS)  =  (NS)  owing  to  the  commutation  of  time  and  spatial  averaging  operations, 
space-time  filtering  and  time-space  filtering  of  the  Navier-Stokes  equation  lead  to  two  dif¬ 
ferent  decompositions  for  the  turbulent  stress  tensor,  a  quantity  that  needs  to  be  modelled 
(viz.,  the  turbulence  closure  problem).  The  subtle  differences  in  these  two  decompositions 
for  the  turbulent  stress  tensor  (arising  from  either  a  space-time  or  time-space  filtering  of  the 
nonlinear  convective  term  in  the  Navier-Stokes  equation)  will  be  elucidated  in  the  following 
subsections. 

Spatial  average  of  the  time-averaged  NS  equation 


The  spatial  average  of  the  time-averaged  NS  equation  (or  spatial  average  of  the  RANS 
equation)  has  been  described  in  detail  by  Raupach  and  Shaw  [14],  Ayotte  et  al.  [25],  and 
others.  Consequently,  only  some  final  results  are  summarized  here  for  later  reference.  The 
spatially-averaged  RANS  equation  for  the  prediction  of  the  spatially-averaged  time-mean 


velocity  (iq)  is 


d(uj)K) 

dxj 


d(P ) 

dxt 


d  .  ,  7 

+  te"(T«)+/" 


(ii) 


1  Note  that  after  spatial  averaging,  a  flow  property  (f)  is  a  continuous  function  of  the  coordinates  in  a  multiply 
connected  space.  Hence,  even  though  <j)  may  vanish  identically  in  the  solid  phase  within  V,  its  spatially-averaged 
value  {<p)e  is  continuous  and  nonzero  in  V,  so  ({<p)e)e  7^  (0)e- 
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with 


and 


Tij  -  -{u'iU'j)  (Ui  Ui)+Udxj’ 

(12) 

f-  =  v!MdS~^Lpn‘dS; 

viscousdrag  formdrag 

(13) 

Here,  p  is  the  kinematic  mean  pressure,  /,  is  the  total  mean  drag  force  per  unit  mass  of 
air  in  the  averaging  volume  composed  of  the  sum  of  a  form  (pressure)  drag  and  a  viscous 
drag,  and  r,j  is  the  spatially- averaged  kinematic  total  stress  tensor.  In  Equation  (13),  v 
is  the  kinematic  viscosity,  S  is  the  part  of  the  bounding  surface  in  the  averaging  volume 
V  that  coincides  with  the  obstacle  surfaces,  rij  is  a  unit  normal  vector  in  the  ith  direction 
pointing  from  V  into  S  (viz.,  directed  from  the  fluid  into  the  solid  surface),  and  dj On  denotes 
differentiation  along  rij.  Note  that  the  spatially-averaged  kinematic  Reynolds  stresses  (u^Uj) 
and  kinematic  dispersive  stresses  (u'/u'j)  are  a  direct  consequence  of  the  spatial  averaging 
of  the  time-averaged  nonlinear  convective  term  UiUy,  viz., 

(uiuj)  =  (ui)(uj)  +  {u'iU'j)  +  {u"u").  (14) 

The  last  term  on  the  right-hand-side  of  Equation  (14)  is  the  dispersive  stress  which  arises 
from  the  spatial  correlation  in  the  time-mean  velocity  field  varying  with  position  in  the 
averaging  volume  V .  Although  Ayotte  et  al.  [25]  rigorously  derive  the  transport  equation 
for  (u'tij),  the  extra  source/sink  terms  in  their  proposed  model  for  the  spatially-averaged 
second  central  velocity  moments  u(u'-,  which  they  denote  as  dij  (contribution  to  the  total 
dissipation  arising  from  the  canopy  interaction  processes),  were  obtained  from  an  approxi¬ 
mate  expression  for  the  work  done  by  the  fluctuating  turbulence  against  the  fluctuating  drag 
force.  The  latter  was  derived  in  the  context  of  the  time  average  of  the  spatially-averaged 
NS  equation.  Mixing  the  spatially-averaged  RANS  formulation  with  the  time-averaged, 
spatially-averaged  NS  formulation  results  in  a  mathematical  inconsistency  in  the  approach 
described  by  Ayotte  et  al.  [25]. 

In  the  next  subsection,  we  formulate  the  equation  set  for  the  time-averaged,  spatially- 
averaged  NS  approach.  The  approach  taken  here  is  similar  to  that  proposed  by  Wang  and 
Takle  [22]  and  Getachew  et  al.  [33]. 

Time  average  of  spatially-averaged  NS  equation 

The  spatial  average  of  the  nonlinear  convective  term  UiUj  in  the  Navier-Stokes  equation  can 
be  expanded  as  follows: 

(utuj)  =  (  «>.,)  +  <)  ((uj)  +  u'J) )  =  («,)<«,)  +  «»">.  (15) 

Using  the  decomposition  for  the  spatially-averaged  nonlinear  convective  term  of  Equa¬ 
tion  (15),  and  applying  the  spatial  averaging  theorem  in  Equation  (2)  with  the  filter  kernel 
defined  in  Equation  (4),  the  spatially-averaged  NS  equation  assumes  the  form 

+  A(t)  +  fu  (i6) 

dxj  oxi  axj 
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where 


and 


Tij  =  —(uiuj)  +  v 


diuj) 

dxj  ' 


(17) 


(18) 


Time-averaging  Equation 
equation  (NS)  as 


(16)  gives  the 

d(uj)(uj)  _  _ 
dxj 


time- averaged ,  spatially-averaged 


d(p) 

dxi 


+ 


d 

dxj 


(Tij )  +  fii 


where 


_  _  _  _  dlui) 

Tij  -  -«)(«')  +  Tij  =  -{«')(«;■)  -  («)  + 


Navier-Stokes 

(19) 

(20) 


and  fi  here  is  the  same  as  fi  defined  in  Equation  (13). 

From  Equations  (11),  (12),  (19)  and  (20),  the  following  relationship  holds: 


+  «*?>  =  K)K-)  +  («)- 


(21) 


Equation  (21)  is  the  necessary  and  sufficient  condition  for  (NS)  =  (NS).  The  total  stress 
tensors  Tij ,  defined  in  either  Equation  (12)  or  Equation  (20)  are  identical  (hence,  the  same 
notation  is  used  for  these  two  quantities),  although  the  individual  terms  in  their  sums  are 
different.  We  note  that  the  physical  character  of  the  term  ( u'-u" )  is  different  from  the 
conventional  dispersive  term  (u'-u"). 


From  the  turbulence  modelling  point  of  view,  we  will  model  (n()(u'  )  in  Equation  (20)  using 
the  Boussinesq  eddy  viscosity  (ut)  closure  as  follows: 


«)K-)  = 

(  d(ui) 
\  dxj 

,  d(uj)\ 
dxi  )' 

(22) 

Here  7t  and  ut  are  defined  as 

«= 

vt  - 

K2 

n  _ 

> 

(23) 

and  e  is  defined  as 

dxk  dxk 

(24) 

In  Equation  (23),  is  a  closure  (empirical)  constant  taken  to  be  0.09  as  in  the  standard 
k-e  model  for  turbulence  closure  [34].  We  note  that  k  is  the  sub-filter  kinetic  energy,  so  k 
is  the  time-averaged  sub-filter  kinetic  energy.  In  addition,  since  the  filter  in  Equation  (5)  is 
positive  (volume  averaging),  k  is  necessarily  a  positive  semi-definite  quantity. 
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To  make  further  progress,  we  assume 


W-Uj)  «  (u'i)(uj) 


(25) 


in  the  present  study  (viz.,  we  will  simply  neglect  the  term  (ii"u")  for  expediency  since 
no  reference  data  exists  at  this  time  to  guide  its  modelling).  Before  we  derive  the  model 
equations  for  7c  and  e,  let  us  examine  the  momentum  sink  ft  (arising  from  the  pressure  and 
viscous  forces  created  by  the  obstacle  elements  in  V)  in  the  spatially-averaged  NS  equation 
[cf.  Equations  (16)  and  (18)].  The  drag  force  term  will  be  parameterized  using  the  following 
common  formulation 

fi  =  -CDA((uj){uj))1/2(ui),  (26) 

where  Cd  is  the  element  drag  coefficient  and  A  is  the  frontal  area  density  (frontal  area  of 
obstacles  exposed  to  the  wind  per  unit  volume).  Equation  (26)  can  be  interpreted  simply 
as  a  definition  for  Cd-  Although  the  drag  coefficient  is  almost  invariably  assumed  to  be  a 
constant  for  a  particular  canopy,  it  will  be  shown  later  that  Cd  is  a  function  of  position 
(, x ,  z)  within  the  canopy. 

Following  from  this  parameterization,  fi  (time-averaged  momentum  sink)  in  Equation  (13) 
is  required  to  be  modelled  (approximated)  as 


fi  =  -CdA^)^))112  (ui). 


(27) 


Substituting  ( rq )  =  (fq)  +  (u[)  into  Equation  (26)  and  using  the  binomial  theorem  to 
approximate  the  square  root  term  (see  also  Getachew  et  al.  [33]),  an  approximate  form  for 
fi  that  is  appropriate  for  time  averaging  can  be  derived  as 


fi  ~  -CdAQ  (  {uj)  H - ^2 - •"  ( ui ) 

(K  )(uj)(u'j)  (Hi  )(u'j)(u'j) 

'  1 


Q2 


2Q2 


(28) 


1  /9 

where  Q  =  ((ui)(ui))  '  is  the  magnitude  of  the  spatially-averaged,  time-mean  wind  speed. 

Time  averaging  of  fi  in  Equation  (28)  gives  the  following  expression  (approximation)  for 
the  time-averaged  form  and  viscous  drag  force  vector  exerted  on  a  unit  mass  of  air  in  the 
averaging  volume: 


h  =  -CoA  (  Q(ui)  + 


(29) 


which,  in  combination  with  Equation  (22),  yields 


h  =  -cda 


r  ,  5  k\,  .  fd(Ui)  d(uj}\  (uj) 

Q+3Q>{Ui)-",{~8fr  ^)~Q 


(30) 
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With  these  closure  assumptions,  the  final  form  of  the  modelled  time-averaged,,  spatially- 
averaged  NS  equation  [obtained  by  substituting  Equations  (22),  (25)  and  (26)  into  Equa¬ 
tions  (19)  and  (20)]  becomes 


d(uj)(uj)  __  _<9(p)  +  _d 


dx-i 


dxi 

CdA 


dxj  L 


( V  +  Vt) 


d(ui)  +  d{uj) 


dxi 


dxi 


2  Sij  K 


^  5  k  \  . 

Q+3 

d{ui)  ,  d(uj)  \  (Uj) 

Q 


-vt 


( +  d^y\ 

V  dxj  dxi  ) 


(31) 


Derivation  of  transport  equations  for  k  and  e 


The  budget  equations  for  k  and  e,  which  have  been  defined  explicitly  in  Equations  (23)  and 
(24),  need  to  be  derived.  To  this  purpose,  let  us  define  f-  =  fa-  fa  as  the  fluctuating  drag 
force.  From  Equations  (28)  and  (29),  we  can  derive 


/'  -  -CdAQ 


'( Uj)(uk)(u'k ) 


Q2 


+  {<)  + 


(n')(^)K) 

Q2 


K) 

Q2 


WfaWk)  + 


(Ui)(u'k)(u'k)  (Uj)K 


2  Q2 


Q2 


(32) 


The  transport  equation  for  the  spatially-averaged  fluctuating  velocity  (u-),  obtained  by 
subtracting  the  evolution  equation  for  the  time-averaged  spatially-averaged  velocity  [Equa¬ 
tion  (19)]  from  the  spatially-averaged  Navier-Stokes  equation  [Equation  (16)],  can  be  written 
in  symbolic  form  as  follows: 

D{v’j)  =  . . .  +  (33) 

Dt  II' 


where  D/Dt  is  the  material  derivative  based  on  the  spatially-averaged  velocity  (ufa.  The 
time  average  of  the  linear  combination  {u'^D^'fa  /  Dt  +  (u') D (u'j) / Dt  gives  the  following 

transport  equation  for  (u£)(u'): 


EMt 

Dt 


+  «> 


D(ur)  =  jK)K-) 

Dt  Dt 


Fiji 


(34) 


where  D/Dt  is  the  material  derivative  based  on  the  spatially-averaged,  time-mean  velocity 
(ui).  Furthermore,  Fy,  representing  the  interaction  between  the  fluctuating  drag  force  and 
spatially-averaged  velocity  fluctuations,  has  the  explicit  form 


Fij  =  (u')fa  +  («•)/• 

2Q{u'i){u'j)  +  ^  ( (ui){uk){u j)(u'k)  +  (uj)(uk){u'i)(u'k)  ) 


=  -CdA 


+  ^(s»>K)(«5)K> 

+^((fii)(«')K)K)  +  <*i)M>K>W>) 


(35) 
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One-half  the  trace  of  Fij  yields 


F  =  -Fa  =  <u')/' 
=  -CdA 


2®k  +  ^  ( (^)(ufc)K)K) ) 


+ 


^((«fc)K)K)K)) 


(36) 


The  triple  correlation  term  (n()  (u()  (u'k)  in  Equation  (36)  can  be  modelled,  following  Daly 
and  Harlow  [35],  as 


K 


K)K)K)  =  2CS^ 


9r 


(u'k)(U'l)Q^i+(U'i)(Ul) 


d«)«) 


dxi 


(37) 


where  the  closure  constant  Cs  ~  0.3  is  used  in  the  present  study.  This  is  a  gradient 
transport  model  for  the  third  moments  of  the  spatially-averaged  fluctuating  velocity  and 
involves  a  tensor  eddy  viscosity.  In  Equations  (36)  and  (37),  the  double  correlation  (Uj)(u'j) 
was  modelled  previously  using  the  constitutive  relationship  in  Equation  (22). 


The  transport  equation  for  the  time- averaged,  sub- filter  kinetic  energy  k  is  obtained  by  mul¬ 
tiplying  Equation  (33)  by  {u[)  and  time  averaging  the  result.  This  procedure  will  give  rise 
to  the  F  term  exhibited  in  Equation  (36).  This  term  represents  the  interaction  of  the  flow 
with  the  obstacle  elements  and  corresponds  explicitly  to  the  work  done  by  the  turbulence 
against  the  fluctuating  drag  force.  The  term  F  can  be  interpreted  as  an  additional  physi¬ 
cal  mechanism  for  the  production/dissipation  of  k  associated  with  work  against  form  and 
viscous  drag  on  the  obstacle  elements.  From  this  perspective,  the  exact  transport  equation 
for  K  is 


where  F  =  (it' )  f- ;  the  flux  Tj  is 

i _  _  pn? 

Tj  =  2  (uj)(ui)(u'i)  +  K)(p'>  - (39) 

and, 

p  =  <4(» 

is  the  production  term  (which  is  generally  positive,  and  hence  a  ‘source-  in  the  k  equation). 
In  addition  to  F,  the  exact  transport  equation  for  k  embodies  an  extra  term  represented 
by  the  second  term  on  the  right- hand-side  of  Equation  (38).  This  term  is  the  energy 
redistribution  due  to  the  interaction  of  the  spatially-averaged  velocity  fluctuations  with  the 
gradient  of  the  sub-filter  stresses  {u"u"). 


The  modelled  transport  equation  for  k  is  then  obtained  as  follows.  Firstly,  the  additional 
energy  redistribution  term  identified  above  is  assumed  to  be  negligible,  and  will  be  ignored 
henceforth.  Secondly,  the  energy  flux  Tj  is  modelled  with  a  gradient  diffusion  hypothesis 


Tj 


vt  9k 

(Tfc  9Xj 


(41) 
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where  the  ‘turbulent  Prandtl  number’  for  k  is  assumed  to  be  ak  =  1.  Thirdly,  the  additional 
physical  effect  on  7c  due  to  viscous  and  form  drag  on  the  obstacle  elements  embodied  in  T1  = 
(u'j )  f  -  will  be  modelled  using  Equations  (36)  and  (37).  With  these  closure  approximations, 
the  model  transport  equation  for  7c  assumes  the  form 


« (^)+(p+F)-e. 

dxj  dxj  \(TkOXjJ 

where  the  explicit  form  for  F  is  exhibited  in  Equations  (36)  and  (37). 


(42) 


The  exact  transport  equation  for  e  can  be  derived  rigorously,  but  it  is  not  a  useful  starting 
point  for  a  model  equation.  Consequently,  rather  than  being  based  on  the  exact  equation, 
the  model  equation  for  e  here  is  essentially  a  dimensionally  consistent  analog  to  the  k- 
equation.  In  this  sense,  the  model  equation  for  e  is  best  viewed  as  being  entirely  empirical. 
To  this  purpose,  we  note  that  the  time  scale  r  =  /c/e  will  make  the  production  and  dissipa¬ 
tion  terms  in  the  7c-equation  dimensionally  consistent.  Hence,  the  dimensionally  consistent 
analog  to  Equation  (42)  becomes 


9{uj)e 

dxj 


+  =(Cei(P  +  F)  -Ct2e), 

Ki 


(43) 


where  at  =  1.3,  Cci  —  1.44  and  Ct 2  =  1.92  are  empirical  (closure)  constants.  The  e-equation 
here  essentially  retains  the  same  form  as  the  usual  model  equation  for  e  commonly  utilized 
in  the  standard  k-e  model  [34].  The  only  difference  here  is  that  the  drag  force  effect  on 
the  turbulence  (embodied  in  the  term  F)  has  been  included  with  the  production  term  P  in 
Equation  (43). 


In  Equation  (43),  we  have  grouped  F  with  P  in  the  e-equation.  In  other  words,  we  sensitize 
the  e-equation  to  the  effects  of  form  and  viscous  drag  of  the  obstacle  elements  by  replacing 
P  with  P  +  F  in  the  ‘production  of  dissipation’  term  (usually,  the  effect  of  the  obstacle 
elements  is  to  enhance  the  dissipation  in  the  canopy  airspace).  This  treatment  is  similar  to 
the  rationale  used  by  Ince  and  Launder  [36]  for  dealing  with  buoyancy  effects  on  turbulence 
in  buoyancy-driven  flows.  In  these  types  of  flows,  the  gravitational  production  term  G  = 
—flgiii'iT'  ( gi  is  the  gravitational  acceleration  vector,  (3  is  the  thermal  expansion  coefficient, 
and  T'  is  the  virtual  temperature  fluctuation)  is  included  with  P  in  the  transport  equation 
for  the  viscous  dissipation  rate.  In  this  regard,  our  proposed  approach  for  treating  F  in  the 
e-equation  differs  from  that  suggested  by  Getachew  et  al.  [33]. 


Wake  production 

The  closure  of  the  spatially- averaged  RANS  equation  [cf.  Equations  (11)  and  (12)]  requires  a 
transport  equation  for  ( k )  =  |(u'u')  (i.e.,  the  spatially-averaged  turbulence  kinetic  energy), 
whereas  that  for  the  time-averaged  spatially-averaged  NS  equation  [cf.  Equations  (19)  and 
(20)]  requires  a  transport  equation  for  k  =  (i.e.,  the  time-averaged  sub-filter 

kinetic  energy).  The  model  transport  equation  for  k  in  Equation  (42)  included  a  source/sink 
term  F  whose  form  can  be  systematically  derived  in  terms  of  the  drag  force  term  that 
appears  in  the  mean  momentum  equation  [cf.  Equation  (36)]. 
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The  budget  equation  for  (k)  can  be  derived  by  applying  the  spatial  averaging  operator  to 
the  standard  transport  equation  for  k  to  give 


d(k )  -r-r\  dui  d 


-(u'jU'iU'i)  +  { u'jP' )  +  -  ( u"  u'X ) 


+  v 


d2(k) 

dx] 


£  —  (  U,U, 


1  3  dxj 


(44) 


II 


In  Equation  (44),  e  =  I/(§yL§yL)  is  the  isotropic  turbulent  dissipation  rate  for  (k). 

The  transport  equation  for  ( k )  contains  two  additional  terms  (designated  I  and  II)  that 
need  to  be  approximated.  Term  I  corresponds  to  the  dispersive  transport  of  ( k )  [analogous 
to  the  dispersive  flux  of  momentum  in  Equations  (11)  and  (12)].  Term  II  can  be  identified 
as  a  wake  production  term  (see  Raupach  and  Shaw  [14])  which  accounts  for  the  conversion 
of  mean  kinetic  energy  to  turbulent  energy  in  the  obstacle  wakes  by  working  of  the  mean 
flow  against  the  drag.  This  term  is  analogous  to  the  F  term  that  appears  in  the  transport 
equation  for  k,  but  unlike  F  whose  form  can  be  systematically  derived  from  the  form 
and  viscous  drag  force  term  that  appears  in  the  mean  momentum  equation,  the  link  (if 
any)  between  the  wake  production  term  in  Equation  (44)  and  the  drag  force  term  in  the 
mean  momentum  equation  is  less  obvious.  For  example,  Raupach  and  Shaw  [14]  showed 
that  provided  (1)  the  dispersive  stress  (u”u”)  and  the  dispersive  transport  of  (k)  are  both 
negligible  and  (2)  the  mean  kinetic  energy  is  not  directly  dissipated  to  heat  in  the  canopy, 
the  wake  production  term  can  be  approximated  as  follows: 

-<K^)"U>  -  =  2 CdAQ  (45) 

3  - 

K 

where  K  represents  the  mean  kinetic  energy  (kinetic  energy  of  the  spatially-averaged  time- 
mean  flow).  Note  that  use  of  Equation  (45)  as  a  model  for  the  wake  production  term  strictly 
provides  a  source  term  in  the  transport  equation  for  ( k ). 


However,  Green  [21]  and  Liu  et  al.  [24]  found  that  it  was  important  to  include  also  a  sink 
term  in  the  budget  equation  for  ( k )  and.  to  this  purpose,  modelled  the  wake  production 
term  in  the  budget  equation  for  ( k )  in  an  ad  hoc  manner  as 

-«“iy)"U>  =  2Cd^Q  QfaXfli))  -4CdAQ  (46) 

3  vj - - - i'  ' - - - " 

K  (fc> 

which  includes  a  gain  to  ( k )  from  conversion  of  the  mean  kinetic  energy  K  to  turbulence 
energy  at  the  larger  scales  (source  term)  and  a  loss  from  ( k )  of  the  large-scale  turbulence  en¬ 
ergy  to  smaller  (wake)  scales  (sink  term).  More  specifically,  Green  [21]  argued  heuristically 
that  the  sink  term  in  Equation  (46)  was  required  to  account  for  the  accelerated  cascade  of 
( k )  from  large  to  small  scales  due  to  the  presence  of  the  roughness  elements  (arising  from 
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the  rapid  dissipation  of  fine-scale  wake  eddies  in  a  plant  canopy).  Liu  et  al.  [24]  noted  that 
the  ad  hoc  inclusion  of  the  sink  term  (second  term  on  the  right-hand  side)  of  Equation  (46) 
[which  was  inserted  in  hindsight]  was  important,  for  otherwise  they  found  that  their  pre¬ 
dicted  (A)  was  “about  100%  larger  than  the  experimental  measurements  when  the  second 
term  was  ignored” . 

In  contrast,  the  additional  source/sink  term  F  that  appears  in  the  budget  equation  for  k  can 
be  systematically  derived  from  rate  of  working  of  the  turbulent  velocity  fluctuations  against 
the  fluctuating  drag  force,  and  appears  naturally  in  the  derivation  of  the  budget  equation 
for  7c.  Even  though  the  ‘turbulence  kinetic  energy’  k  =  (time- averaged,  sub-filter 

kinetic  energy)  used  in  our  turbulence  closure  model  is  different  from  the  usual  form  of  the 
spatially-averaged  turbulence  kinetic  energy  (A)  =  \{ u [u^),  they  are  nevertheless  related 
as  follows:  1  _ _ 

<*>  -R  =  !(<<<)  -  <«> )  =  2<  K)"«)"  )•  (47) 

Note  that  the  difference  between  (A)  and  7c  is  proportional  to  the  difference  between  the 
two  forms  of  dispersive  stress  that  appear  in  the  spatially-averaged  RANS  equation  and 
the  time-averaged,  spatially-averaged  NS  equation.  However,  note  that  this  difference  can 
be  expressed  as  the  spatial  average  of  time  averages  of  the  departures  of  velocity  fluctu¬ 
ations  from  their  spatial  (volume)  average.  Since  this  term  involves  a  “perturbation  of  a 
perturbation” ,  it  seems  reasonable  to  assume  that 

0«  (K)"K)" )  |  «max((A),7c).  (48) 

A 

With  this  assumption,  (A)  and  7c  are  expected  to  be  almost  equal  in  value  (viz.,  (A)  «  k). 
This,  together  with  Equation  (21),  also  implies  that 

<u"u")  «  <<<>,  (49) 

or,  in  other  words,  the  dispersive  stresses  are  expected  to  be  approximately  equal  to  the 
spatial  average  of  the  high-frequency  turbulent  stresses.  Finally,  with  reference  to  Equa¬ 
tion  (21),  this  implies  that  (ii')(w')  ~  («(«'• ).  The  latter  approximation  will  be  used  in  a 

later  section  to  compare  model  predictions  of  {u[){u'j)  with  the  diagnosed  values  of  ( U(Uj ) 
obtained  from  a  high-resolution  RANS  simulation. 

Interestingly,  the  ‘zeroth-order’  term  in  our  expansion  of  F  in  Equation  (36)  rewritten  below 
F  =  Wf[  =  —2 CdAQ  Q  («')<«'))  +H.O.T, 


where  H.O.T  denotes  higher-order  correction  terms,  is  analogous  to  the  sink  term  contribu¬ 
tion  for  the  wake  production  term  in  Equation  (46)  [except  the  factor  here  is  2,  rather  than 
4],  Note  that  the  higher-order  correction  terms  for  F  can  have  either  sign  implying  that  they 
are  source/sink  terms.  More  importantly,  it  needs  to  be  emphasized  that  the  leading-order 
term  of  F  is  a  sink  term,  and  not  a  source  term  implying  that  in  the  transport  equation  for 
k  the  conversion  of  MKE  to  TKE  is  ipso  facto  absent. 
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Determination  of  drag  coefficient 


In  this  section,  we  will  explicitly  diagnose  a  drag  coefficient  Co  using  the  results  for  a 
high-resolution  CFD  simulation  of  a  developing  flow  over  an  aligned  array  of  cubes  (3-D 
buildings)  (described  in  I).  This  regular  building  array  consisted  of  seven  rows  of  cubes 
arranged  with  a  plan  area  density  of  0.25  [28].  To  derive  a  drag  coefficient  for  this  array, 
first  note  from  Equations  (13)  and  (29)  that  the  ^-component  of  the  drag  force  fa  (i  =  1  =  x) 
reduces  to 

/.  s  £  jf  -  ^JsPn,dS  =  -CD(z)AQ(u),  (51) 

where  we  assume  implicitly  that  V  in  Equation  (51)  is  a  thin  horizontal  slab  centered  at 
the  level  z. 

To  specialize  the  analysis  to  the  aligned  3-D  building  array  consisting  of  the  7  rows  of  cubes 
arranged  normal  to  the  mean  wind  direction,  we  consider  decomposing  this  array  into  seven 
“drag  units"  as  illustrated  in  Figure  1(a).  Each  drag  unit  consists  of  one  row  of  cubes 
(buildings)  and  the  associated  downstream  street  canyon.  To  diagnose  a  drag  coefficient 
Cp  (z)  at  level  z  for  each  drag  unit  in  the  array,  we  take  an  averaging  volume  V  as  the  thin 
horizontal  slab  shown  in  Figure  1(b).  Now  application  of  Equation  (51)  to  this  averaging 
volume  gives  the  following  explicit  expression  for  the  pressure  drag: 

If  1 

—  y  J  pnxdS  =  ~2jp  J  {p\x=xO  ~  P\x=xO+h)  dy  .  (52) 

The  frictional  drag,  y  fs  §%dS,  can  be  evaluated  similarly  using  the  ‘wall  function’  approach 
to  estimate  the  skin  friction  along  the  building  surfaces  [37].  With  the  computation  of  the 
pressure  and  frictional  drag  forces,  the  local  or  sectional  drag  coefficient  Cx>(z)  can  then  be 
diagnosed  as  follows  [cf.  Equation  (51)]: 


Cd{z)A 


-fx 

max(  5,  Q(u)(z) )  ’ 


(53) 


where  6  ss  0.0025  <C  1  is  chosen  here  to  avoid  a  possible  singularity  problem  for  the  case 
when  ( u )  «  0.  We  note  that  the  results  to  be  shown  in  the  next  section  are  insensitive 
to  the  exact  value  chosen  for  6  in  Equation  (53).  Furthermore,  for  developing  flow  over  a 
building  array,  a  local  equilibrium  assumption  involving  the  balance  between  the  obstacle 
drag  force  and  the  downward  transport  of  momentum  by  turbulent  stresses  cannot  be  used 
to  diagnose  the  sectional  drag  coefficient. 


Vertical  profiles  of  Cx>(z)  over  the  range  of  heights  0  <  z/H  <  1  for  each  drag  unit  in 
the  3-D  building  array  are  shown  in  Figure  2.  For  the  first  drag  unit  (associated  with  the 
first  row  of  buildings  in  the  array),  Cd(z)  increases  monotonically  with  decreasing  z/H. 
For  drag  units  in  the  array  interior  (drag  unit  #3  or  greater),  Cd(z)  first  increases,  then 
decreases,  and  then  increases  again  with  decreasing  z/H.  In  all  cases,  the  highest  values  of 
Cd{z )  were  obtained  near  the  ground  as  z/H  — >  0,  with  Cx>(z )  decreasing  to  a  small  value 
at  the  top  of  the  canopy,  the  latter  of  which  implies  a  significant  reduction  in  the  rearward 
pressure  deficit  near  the  canopy  top. 
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We  note  that  Cd  is  smaller  for  the  drag  units  in  the  interior  of  the  building  array  compared 
to  that  of  the  first  drag  unit  which  is  influenced  by  the  impingement  region  upwind  of  the 
array.  This  is  an  example  of  the  ‘shelter  effect’  where  the  mutual  interference  of  the  obstacle 
elements  leads  to  a  reduction  in  the  value  of  the  drag  coefficient  measured  in  situ  within 
the  array  compared  to  the  drag  coefficient  measured  in  the  free-stream  flow  (e.g.,  drag 
coefficient  corresponding  to  drag  unit  #1)-  Indeed,  the  largest  value  of  the  sectional  drag 
coefficient  Cp(z)  at  a  given  z  occurs  in  drag  unit  #1  (i.e..  first  row  of  buildings)  over  much 
of  the  canopy  depth.  This  first  drag  unit  is  responsible  primarily  for  the  adverse  pressure 
gradient  in  the  impingement  zone  upwind  of  the  building  array  that  results  in  a  significant 
deceleration  of  the  air  flow  here.  In  contrast,  the  smallest  value  of  the  drag  coefficient 
Cx>(z)  at  each  level  z/H  occurs  for  drag  unit  #2  which  corresponds  to  the  region  in  the 
building  array  (viz.,  in  the  immediate  leeward  zone  of  the  first  row  of  buildings)  where  the 
shelter  effect  is  most  pronounced.  More  specifically,  the  maximum  mean  wind  speed  and 
turbulence  reduction  occurs  in  this  immediate  leeward  ‘quiet’  zone  (corresponding  to  drag 
unit  #2). 


Because  such  detailed  information  on  the  sectional  (local)  drag  coefficient  is  not  usually 
available,  it  is  useful  to  consider  also  the  application  of  a  bulk  drag  coefficient.  To  this 
purpose,  we  consider  three  different  methods  for  the  evaluation  of  a  bulk  drag  coefficient, 
Cd  bulk-  These  three  methods  for  determination  of  a  bulk  drag  coefficient  are  defined  as 
follows: 

Method  1:  z_H 

Cd, bulk, i  =  T7  f  Co(z)dz ,  (54) 

zi  Jz= o 

where  Cd, bulk/  (1  <  f  <  7)  is  the  bulk  drag  coefficient  for  each  drag  unit  determined  by 
averaging  the  sectional  drag  coefficient  at  height  z  (Cd(z))  over  the  entire  depth  of  the 
canopy; 


Method  2: 


1  7 

Cd, bulk,ave  =  j  ^  ^  Cp.bulk.n 

i=  1 


(55) 


is  the  arithmetic  average  of  Cd, bulk/  over  the  drag  units  comprising  the  entire  3  D  building 
array;  and, 


Method  3: 

CD,bulk,l  =  CD,bulk/|i=l 


(56) 


is  the  bulk  drag  coefficient  for  the  first  drag  unit.  This  is  similar  to  using  a  drag  coefficient 
diagnosed  for  the  flow  over  a  single  (isolated)  row  of  buildings  (viz.,  a  free-stream  bulk  drag 
coefficient). 


The  variation  of  the  bulk  drag  coefficient  Cd, bulk/  (defined  using  Method  1  above)  with 
increasing  drag  unit  i  in  the  array  (viz.,  with  increasing  fetch  through  the  urban  canopy) 
is  exhibited  in  Figure  3(a).  Here,  we  have  decomposed  the  drag  coefficient  Cd, bulk/  Into 
a  contribution  from  the  form  (pressure)  drag  Cdp, bulk/  and  the  frictional  (viscous)  drag 
Cd/, bulk/-  Note  that  the  contribution  of  the  frictional  (viscous)  drag  (Cd/, bulk)  to  the 
bulk  drag  coefficient  is  negligible  in  comparison  with  the  contribution  from  the  pressure 
(form)  drag  Cdp, bulk-  Generally,  the  contribution  of  the  skin  friction  to  the  bulk  drag 
coefficient  is  less  than  one  percent  of  the  contribution  of  the  form  drag.  We  note  that  the 
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bulk  drag  coefficient  for  drag  units  in  the  array  interior  (viz.,  for  drag  unit  #3  or  greater) 
is  between  one-half  and  two-thirds  of  the  bulk  drag  coefficient  for  drag  unit  #1  (viz.,  drag 
unit  associated  with  the  row  of  buildings  forming  the  upstream  boundary  of  the  array) 
implying  a  ‘shelter  factor’  (ratio  of  free-stream  bulk  drag  coefficient  to  the  in  situ  value) 
of  between  about  1.5  and  2.  This  is  consistent  with  the  reported  range  of  shelter  factors 
of  between  1  to  4  compiled  by  Raupach  and  Thom  [38]  in  their  review  of  canopy  flows. 
Finally,  the  bulk  drag  parameters  Cd, bulk, aveA  and  Cx>,buik,iA  were  found  to  be  about  0.7 
and  1.2,  respectively  [cf.  Figure  3(b)]  for  the  3-D  building  array  considered  here. 

Results  and  discussion 


The  results  to  be  presented  in  this  section  were  obtained  with  three  different  versions  of  the 
k-e  model  (distinguished  by  the  source/sink  terms  included  in  the  transport  equations  for 
k  and  e)  with  the  major  differences  between  these  versions  summarized  below.2 

Model  (Version)  1:  Model  1  is  in  essence  the  same  as  the  proposed  model  except  that  /*  in 
Equation  (30)  is  replaced  by 

fi  =  -CDAQ{v,i),  (57) 

and  F  in  Equation  (36)  is  set  identically  to  zero.  Hence,  in  this  version  of  the  model,  no 
additional  source  or  sink  terms  are  included  in  the  transport  equations  for  k  and  e,  and  a 
simplified  parameterization  for  the  mean  drag  force  vector  [viz.  Equation  (57)]  is  used. 

Model  (Version)  2:  Model  2  is  Lee  and  Howell’s  model  [39],  in  which  the  parameterization 
for  the  mean  drag  force  vector  fi  is  that  displayed  in  Equation  (57).  However,  unlike 
Model  1,  the  leading  order  term  of  F  in  Equation  (36)  is  retained  (viz.,  a  sink  term  of 
the  form  F  «  —2CdAQk  is  incorporated  into  the  transport  equation  for  k).  Finally, 
rather  than  include  this  leading  order  approximation  of  F  with  the  production  term  P  as  in 
Equation  (42),  Lee  and  Howell’s  model  introduces  a  sink  term  of  the  form  eF/R  ~  —2CpAQe 
into  the  transport  equation  for  e.  However,  we  note  that  the  inclusion  of  the  ‘zeroth- 
order’  approximation  for  F  here  with  the  production  term  P  in  the  e-transport  equation 
[cf.  Equation  (42)]  did  not  make  a  significant  difference  in  the  flow  predictions  (viz.,  the 
predictions  of  urban  flow  using  the  latter  model  are  indistinguishable  from  those  obtained 
using  Lee  and  Howell’s  model). 

Model  (Version)  3:  Model  3  is  the  current  model  (described  above),  consisting  of  Equa¬ 
tions  (31),  (36),  (37),  (40),  (42),  and  (43). 


The  K-isopleths  predicted  by  three  different  models  (summarized  above)  using  the  local 
(sectional)  drag  coefficient  Cp{z)  (cf.  Figure  2)  are  displayed  in  Figure  4.  As  the  upstream 
flow  impinges  on  the  leading  edge  of  the  building  array,  which  is  modelled  here  as  a  ‘porous 
layer’  of  height  H  extending  from  the  upstream  edge  of  the  first  drag  unit  to  the  downstream 
edge  of  the  last  drag  unit,  the  deceleration  of  the  flow  within  the  volume  of  the  canopy 
{z/H  <  1)  and  acceleration  (overspeeding)  of  the  flow  above  about  z/H  k.  4/3  results  in  a 

2  Although  we  use  the  standard  terminology  in  k-e  closure  modelling  and  refer  to  the  transport  equation  for  the 
turbulence  energy  as  the  fc-equation.  this  should  be  interpreted  more  precisely  here  to  refer  actually  to  the  transport 
equation  for  K  (viz.,  the  time-averaged,  sub-filter  kinetic  energy).  Consequently,  in  this  section,  the  ^-equation 
will  be  interpreted  to  be  synonymous  with  the  K-equation,  so  that  k  and  K  will  be  used  interchangeably  in  what 
follows. 
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strong  shear  in  the  mean  flow  in  the  region  of  the  canopy  top.  An  inflection  point  occurs 
in  the  mean  wind  at  z/H  ~  1  where  the  mean  shear  d(u)/dz  assumes  its  maximum  value. 
This  mean  shear  acts  as  a  source  for  the  generation  of  k  through  P  ~  ( d(u)/dz )  (in  the 
Boussinesq  eddy-viscositv  approximation)  at  or  near  the  top  of  the  canopy.  The  presence  of 
an  inflection  point  in  the  mean  velocity  profile  and  the  strong  generation  of  k  in  the  region 
of  the  inflection  point  is  characteristic  of  a  plane  mixing  layer  (i.e.,  a  free  shear  layer  that 
forms  when  two  airstreams  of  different  velocities,  initially  separated,  merge  downstream 
of  the  separation  and  interact)  as  discussed  by  Raupach  et  al.  [40]  for  the  case  of  plant 
canopies. 


A  perusal  of  the  results  in  Figure  4  shows  that  Model  1  performs  significantly  differently  than 
the  other  two  models.  In  particular,  Model  1  exhibits  a  much  more  prominent  turbulence 
kinetic  energy  (TKE)  ‘plume’  than  is  evident  in  the  other  2  models.  This  TKE  ‘plume’ 
originates  from  the  upstream  rooftop  edge  of  the  first  building  and  spreads  vertically  at  a 
constant  (approximately  or  better)  rate,  achieving  a  vertical  depth  of  about  5 H  at  roughly 
a  downstream  distance  of  25 H  from  the  upstream  edge  of  the  building  array.  The  TKE 
‘plume’  is  advected  by  the  fluid  flow  in  the  downstream  direction,  and  its  rate  of  spread  in  the 
vertical  direction  is  governed  by  turbulence  diffusion  through  the  eddy  viscosity  v%  k  .  For 
Model  1,  no  drag- force-related  ‘sink’  terms  are  included  in  the  k-  and  e-equations,  resulting 
in  an  excessive  level  of  k  and.  hence,  Vf  This,  in  turn,  increases  the  spreading  rate  of  the 
TKE  ‘plume’.  In  contrast,  the  turbulence  energy  levels  predicted  by  Models  2  and  3  (which 
incorporate  minimally  an  additional  sink  term  in  the  k-  and  e-equations)  are  smaller  than 
those  predicted  by  Model  1.  Furthermore,  Model  1  also  predicts  a  much  higher  level  of 
k  within  the  volume  of  the  canopy  (particularly  in  the  near  wall  region  0  <  z/H  <  0.5) 
than  Models  2  or  3.  In  general,  the  prediction  of  turbulence  energy  provided  by  Models 
2  and  3  is  similar  except  for  the  region  downstream  of  the  leeward  rooftop  edge  of  the 
last  building  (14  <x/H<  17.5)  where  Model  3  predicts  a  slightly  higher  level  of  k  than 
Model  2.  It  appears  from  this  observation  that  k  is  not  very  sensitive  to  the  differing 
treatments  of  sources  and  sinks  in  the  k-  and  e-equations,  while  simulations  using  the  k-e 
closure  without  additional  sources  and  sinks  in  the  k-  and  e-equations  lead  to  significantly 
different  predictions  of  the  turbulence  energy  levels. 


Profiles  of  predicted  spatially-averaged  time-mean  streamwise  velocity  ( u )  obtained  using 
Models  1-3  for  the  first  six  drag  units  are  shown  in  Figure  5.  These  predictions  are  compared 
to  (u)  obtained  by  spatially-averaging  RANS  mean  flow  solutions  over  the  same  volumes 
occupied  by  the  drag  units.  Our  high-resolution  RANS  simulations  of  the  3-D  building  array 
considered  here  are  reported  in  I.  The  latter  results,  which  have  been  spatially-averaged 
here,  will  be  labelled  in  the  figures  using  the  symbol  ‘o’  and  will  be  referred  to  henceforth 
as  ‘high-resolution  CFD’  results.  Within  the  volume  of  the  canopy,  the  spatially-averaged 
time-mean  velocity  predicted  by  Models  2  and  3  (cf.  Figure  5)  are  in  excellent  agreement 
with  the  high-resolution  CFD  results  for  drag  units  #3  to  #6.  However,  for  the  first  two 
drag  units,  (u)  is  slightly  overpredicted  by  these  two  models  which  appears  to  be  a  result 
of  the  severe  underestimation  of  (u')(u')  in  the  near-wall  region  z/H  <  2/3  (see  below). 
Vertical  profiles  of  ( u )  obtained  from  the  high-resolution  CFD  results  appear  to  have  reached 
streamwise  equilibrium  (viz.,  the  spatially-averaged,  time-mean  streamwise  velocity  field  is 
fully  developed)  by  the  third  drag  unit,  a  feature  that  is  reproduced  correctly  by  all  three 
models. 
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The  corresponding  profiles  of  turbulence  quantities  ( u')(u ')  (streamwise  normal  stress)  and 
( u'){w ')  (shear  stress)  are  displayed  in  Figures  6  and  7,  respectively.  Consistent  with  Fig¬ 
ure  4  in  relation  to  the  spreading  of  the  TKE  ‘plume’,  the  peak  values  of  (u')(u')  (i.e. ,  the 
streamwise  normal  stress  component  of  k )  at  z/H  ft  4/3  are  overpredicted  by  all  three 
models.  This  overprediction  of  ( u'){u ')  is  rather  significant  for  Model  1,  particularly  in 
drag  units  #2  to  #6.  An  overestimation  of  { u ,){u')  (and,  by  inference  of  k )  gives  too  high 
values  of  the  eddy  viscosity  vt  ~  k 2,  yielding  an  excessive  vertical  turbulent  diffusion  of 
mean  momentum.  This,  in  turn,  results  in  an  underestimation  of  the  magnitude  of  the 
mean  shear  at  the  top  of  the  canopy  and  an  underprediction  of  ( u )  in  the  region  above  the 
canopy  height  (i.e.,  in  z/H  >  1). 

The  vertical  spreading  of  the  TKE  ‘plume’  is  evident  on  an  examination  of  the  streamwise 
development  of  ( u')(u ')  in  Figure  6.  Here,  it  is  seen  that  (u')(u')  approaches  its  undisturbed 
far  upstream  value  at  the  level  of  z  &  350  mm  {z/H  ft  2.3)  in  drag  unit  #1,  but  at  the 
greater  level  of  z  a  500  mm  ( z/H  ft  3.3)  in  drag  unit  #6.  Interestingly,  the  peak  value  of 
(u')(u')  also  increases  monotonically  in  the  streamwise  direction,  which  is  markedly  different 
from  what  we  normally  would  expect  from  the  development  of  a  mixing  layer  (assuming 
that  the  shear  flow  near  the  canopy  top  has  the  characteristics  of  a  plane  mixing  layer). 
For  example,  in  Model  1,  the  peak  value  of  {u')(u')  is  ft  0.3  m2  s“2  in  drag  unit  #1  and 
ft  0.5  m2  s~2  in  drag  unit  #6.  However,  the  high-resolution  CFD  results  show  almost  no 
tendency  of  streamwise  development  in  the  streamwise  normal  stress  after  drag  unit  #3, 
which  is  consistent  with  the  streamwise  variation  of  C/ybuik  displayed  previously  in  Figure  3. 
In  particular,  within  the  array  interior  {x/H  >  6),  it  appears  that  the  flow  statistics  (mean 
flow  and  turbulence)  within  the  third  street  canyon  should  be  “typical”  of  their  neighbors, 
implying  that  in  the  array  interior  a  periodic  boundary  condition  could  have  been  imposed 
to  model  the  flow  within  a  unit  drag  cell  consisting  of  a  single  row  of  3-D  buildings  and  the 
associated  street  canyon. 

The  deficiency  of  not  being  able  to  predict  a  fully-developed  state  for  the  turbulence  fields 
within  the  array  using  a  distributed  drag  force  approach  might  be  due  to  the  defects  in  the  F 
term  included  in  the  transport  equation  for  k.  Although  in  the  present  model,  higher-order 
corrections  to  the  nonlinear  drag  force  term,  similar  to  those  proposed  by  Getachew  et  al. 
[33]  have  been  included,  their  overall  effect  on  the  model  prediction  does  not  appear  to  be 
significant.  Indeed,  the  model  predictions  obtained  with  the  more  complex  model  proposed 
here  is  not  significantly  different  from  the  simpler  Lee  and  Howell  model  [39].  However,  the 
effectiveness  of  the  higher-order  corrections  in  our  model,  particularly  the  triple  correlations 
(u1/)  {u[) (u'k) ,  can  only  probably  be  realized  if  the  Reynolds  stress  anisotropy  can  be  correctly 
represented,  which  is  difficult  to  achieve  within  the  present  impoverished  k-e  modelling 
framework  used  in  conjunction  with  the  Boussinesq  linear  stress-strain  relationship. 


As  the  major  component  in  P  is  —{u,)(w')d(u) /dz,  we  expect  the  shear  stress  —(u'){w')  to 
be  overestimated  by  the  three  models  in  order  to  be  consistent  with  the  overestimation  of  k 
and  (u')(u')  observed  in  Figures  4  and  6.  It  is  seen  in  Figure  7  that  Model  1  predicts  a  peak 
value  in  shear  stress  that  is  approximately  4-5  times  smaller  (i.e.,  more  negative)  than  the 
high-resolution  CFD  results  for  drag  units  #1  to  #6.  Even  for  Models  2  and  3,  the  peak 
values  of  ( u')(w ')  are  still  too  small  (i.e.,  too  negative),  although  Model  2  performs  slightly 
better  than  Model  3  in  the  prediction  of  ( u'){w ')  and  {u'){u')  (cf.  Figure  6).  However,  it 
should  be  emphasized  here  that  the  true  performance  of  the  present  mode!  can  only  be 
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assessed  if  a  higher-order  closure  model,  such  as  the  Reynolds  stress  transport  model  of 
Launder  et  al.  [41]  that  is  capable  of  predicting  the  anisotropy  of  the  turbulence  fields,  is 
used.  This  is  an  interesting  topic  for  a  future  investigation. 

The  effect  of  the  specification  of  Co  on  the  predictive  performance  of  a  distributed  drag 
force  model  will  be  investigated.  To  this  purpose,  the  diagnosed  sectional  drag  coefficient 
(cf.  Figure  2)  and  three  bulk  drag  coefficients  (cf.  Figure  3)  have  been  used  as  input  to 
Model  3  for  the  prediction  of  (u),  { u'){u' ),  and  (u')(w').  The  results  of  this  simulation 
are  displayed  in  Figures  8,  9,  and  10.  It  is  seen  from  Figure  8  that  the  use  of  a  sectional 
drag  coefficient  Cd(z )  permitted  a  more  accurate  prediction  of  (u)  in  drag  unit  #1  for  the 
region  z/H  <  2/3.  This  is  the  result  of  the  correct  representation  in  the  increase  of  the 
drag  coefficient  with  decreasing  zlH  which  is  embodied  in  the  local  drag  coefficient,  but 
is  missing  in  the  representation  provided  by  the  bulk  drag  coefficients.  Nevertheless,  the 
predictive  performance  of  the  model  using  only  the  bulk  drag  coefficients  generally  gave 
similar  results  to  those  obtained  with  the  sectional  drag  coefficient  for  both  the  mean  flow 
and  turbulence. 

For  drag  unit  #2,  it  is  observed  that  the  model  prediction  for  ( u )  below  the  canopy  height 
obtained  using  the  bulk  drag  coefficient  Co,b uik,i  is  markedly  different  from  the  other  three 
sets  of  predictions  obtained  using  Cp(z),  Cd, buik,i>  and  Cr>,buik,ave-  More  specifically,  (u) 
is  too  small  (i.e.,  underpredicted)  for  2/5  <z/H<  1  compared  to  the  high-resolution  CFD 
results.  This  is  because  C^buiM  can  be  regarded  as  the  drag  coefficient  for  a  flow  over 
an  isolated  row  of  cubical  buildings  (i.e.,  a  free-stream  drag  coefficient).  The  drag  of  an 
isolated  row  of  buildings  is  different  from  its  effective  drag  in  situ  (viz.,  the  drag  of  a  row 
of  buildings  within  a  building  array,  where  the  ‘shelter  effect’  becomes  important).  It  was 
shown  earlier  that  the  drag  for  a  row  of  buildings  (modelled  as  a  porous  barrier  here)  in 
situ  is  less  than  for  the  same  row  of  buildings  in  isolation  owing  to  the  shelter  effect  .  In 
fact,  the  current  calculations  show  that  the  drag  coefficient  drops  considerably  after  the 
first  row  of  buildings  in  the  array  due  to  the  ‘shelter  effect’  and  recovers  rapidly  thereafter, 
reaching  approximately  its  fully-developed  value  somewhere  between  the  third  and  fourth 
rows  of  buildings  in  the  array  (cf.  Figure  3). 


Overall,  good  predictions  of  the  mean  flow  through  the  canopy,  except  in  drag  unit  #1, 
can  be  achieved  by  the  using  Cx>,buik,avo  the  latter  of  which  can  be  considered  to  be  equal 
(approximately  or  better)  to  the  fully-developed  value  of  Cx>jbuik  *n  the  array  interior.  Since 
CD(z)  and  CD,buik,i  are  difficult  to  obtained  in  practice,  we  recommend  the  approach  of 
using  C/ybuik.i  (free-stream  drag  coefficient)  for  drag  unit  #1  and  Cp;buik,ave  f°r  all  other 
drag  units  in  the  array  (all  of  which  will  ‘experience’  some  form  of  sheltering  from  the  drag 
unit(s)  upstream  of  it).  Both  Cx>ibuik,i  and  C/ybuik.ave  are  simpler  to  obtain  experimentally 
than  the  sectional  drag  coefficient  Cd{z)- 

Conclusions  _ 


The  present  k-e  model  was  derived  rigorously  based  on  a  time-averaging  of  the  spatially- 
averaged  NS  equation.  If  we  assume  that  the  kinematic  dispersive  stresses.  (u"n"),  are 
comparable  in  magnitude  to  (u"u")  [an  explicit  demonstration  of  the  validity  of  this  as¬ 
sumption  will  require  either  a  Direct  Numerical  Simulation  (DNS)  or  a  Large  Eddy  Simula¬ 
tion  (LES)],  then  (u')(u')  can  be  compared  directly  with  the  spatially-averaged  kinematic 
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Reynolds  stresses,  (u'u)),  the  latter  of  which  can  be  estimated  explicitly  using  the  high- 
resolution  CFD  results  presented  in  I  (or,  equivalently,  7c  can  be  compared  directly  to  (k), 
the  latter  quantity  of  which  is  available  directly  from  a  RANS  simulation). 


We  found  that  the  peak  of  the  kinematic  streamwise  normal  stress,  (u'){u'),  was  overpre¬ 
dicted  by  all  three  models  investigated  here.  However,  we  note  that  these  models  overpredict 
( u')(u ')  (and,  by  inference,  7c)  when  compared  to  the  high-resolution  CFD  results,  the  latter 
of  which  are  assumed  implicitly  to  correspond  to  the  “ground  truth” .  However,  this  may  not 
necessarily  be  the  “correct”  comparison  since  it  was  shown  in  I  that  the  RANS  simulation 
generally  gave  an  underprediction  of  the  turbulence  energy  in  the  building  array,  implying 
that  the  spatially-averaged  TKE  obtained  from  these  high-resolution  CFD  simulations  will 
probably  underpredict  the  true  (but  unknown)  values.  Hence,  the  overprediction  of  the 
streamwise  normal  stress  by  the  distributed  drag  force  models  may  not  be  as  severe  as 
apparently  indicated  here  in  comparison  with  the  high-resolution  CFD  results. 

Nevertheless,  Model  1  (i.e.,  model  with  no  additional  source/sink  terms  in  the  k-  and  e- 
equations)  gave  the  worst  predictive  performance  of  the  three  models  considered.  Model  2 
(Lee  and  Howell’s  model)  and  Model  3  (current  proposed  model)  gave  very  similar  ^predic¬ 
tions  for  the  mean  streamwise  velocity  (u).  In  terms  of  turbulence  quantities  (viz.,  (u'}{u') 
and  ( u')(w ')),  Model  2  gave  slightly  better  results  than  Model  3,  which  appears  to  suggest 
that  the  contribution  of  the  higher-order  correction  terms  in  the  present  model,  displayed 
in  the  F  term  of  Equation  (36),  is  not  important.  In  other  words,  the  predictions  of  the 
turbulence  quantities  do  not  appear  to  be  sensitive  to  the  precise  treatment  of  the  source 
and  sink  terms  in  the  k-  and  e-equations.  However,  it  is  possible  that  the  insensitivity  of 
the  results  to  the  inclusion  of  such  higher-order  terms  in  the  model  is  caused  by  the  use  of 
the  Boussinesq  linear  stress-strain  constitutive  relation  [i.e.,  Equation  (22)]  within  the  k-e 
modelling  framework.  Hence,  higher-order  turbulence  closure  models,  such  as  the  Reynolds 
stress  transport  model  of  Launder  et  al.  [41]  that  is  capable  of  predicting  the  anisotropy 
in  the  Reynolds  stress  field,  might  be  required  here  in  order  to  take  full  advantage  of  the 
higher-order  correction  terms. 

The  profiles  of  ( u )  are  in  general  well  predicted  by  the  present  model  in  conjunction  with 
the  use  of  a  sectional  drag  coefficient  Cd{z ),  or  the  bulk  drag  coefficients  Co.buik.;  and 
C'd, bulk. avej  particularly,  below  the  canopy  height  ( z/H  <  1).  However,  when  Cx^bulk.i  was 
used,  the  (u) -profiles  within  the  canopy  were  generally  underpredicted  (at  least  for  the  first 
four  drag  units  of  the  array)  due  to  the  neglect  of  the  ‘shelter  effect’  in  the  canopy  interior. 
In  consequence,  we  recommend  that  in  modelling  of  the  developing  flow  over  a  building 
array  using  the  distributed  drag  force  approach  that  C^buiM  (or,  equivalently,  the  free- 
stream  drag  coefficient)  be  assigned  to  the  first  drag  unit  (e.g.,  first  row  of  buildings  in  the 
array)  and  Cx>,buik,ave  (or,  equivalently,  an  in  situ  drag  coefficient  for  the  fully-developed 
flow  in  the  array  interior)  be  assigned  to  the  remaining  drag  units. 
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Figure  1.  Decomposition  of  the  aligned  array  of  cubes  into  (a)  seven  drag  units  with  each 
drag  unit  consisting  of  a  row  of  cubes  and  the  associated  downstream  canyon;  (b)  thin 
horizontal  slab  chosen  as  the  averaging  volume  within  each  drag  unit. 
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Figure  2.  The  variation  of  normalized  sectional  (local)  drag  coefficient  Cd(z)AH/V 
Cd(z)A  for  each  drag  unit  comprising  the  building  array. 
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(a) 


(b) 


Figure  3.  The  normalized  bulk  drag  coefficients  for  various  drag  units  in  the  build¬ 
ing  array:  (a)  Cd, bulk, %AH/V ,  C  Dp, bulk, iAH  /V  and  C d  f  ,bu\k,iAH  /V ,  (b)  Cd, bulk, iAH/V , 

C D, bulk, avc  AH /V ,  and  C D, bulk, \ AH /V . 
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Figure  4.  Turbulence  kinetic  energy  ( k )  isopleths  obtained  with  three  different  k-e  closure 
models  used  in  conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient 

CD(z). 
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Figure  5.  Profiles  of  (u)  obtained  with  three  different  k-e  closure  models  used  in  conjunc¬ 
tion  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  Cd(z)  for  drag  units 
#1  and  #2.  Model  1  (solid  line);  Model  2  (dashed  line);  Model  3  (dashed-dotted  line); 
high-resolution  CFD  (open  circle). 
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Figure  5.  (Continued)  Profiles  of  (u)  obtained  with  three  different  k-e  closure  models  used 
in  conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  Cp(z)  for 
drag  units  #3  and  #4.  Model  1  (solid  line);  Model  2  (dashed  line);  Model  3  (dashed-dotted 
line);  high-resolution  CFD  (open  circle). 
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Figure  5.  (Continued)  Profiles  of  ( u )  obtained  with  three  different  k-e  closure  models  used 
in  conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  Cp(z)  for 
drag  units  #5  and  #6.  Model  1  (solid  line);  Model  2  (dashed  line);  Model  3  (dashed-dotted 
line);  high-resolution  CFD  (open  circle). 
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Figure  6.  Profiles  of  (u/)(u,)  obtained  with  three  different  k-c  closure  models  used  in 
conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  Cd(z)  for 
drag  units  #1  and  #2.  Model  1  (solid  line);  Model  2  (dashed  line);  Model  3  (dashed-dotted 
line);  high- resolution  CFD  (open  circle). 
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Unit  #3 


Figure  6.  (Continued)  Profiles  of  ( u'){u ')  obtained  with  three  different  k-e  closure  models 
used  in  conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  C d  (z) 
for  drag  units  #3  and  # 4.  Model  1  (solid  line);  Model  2  (dashed  line);  Model  3  (dashed- 
dotted  line);  high-resolution  CFD  (open  circle). 
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Figure  6.  (Continued)  Profiles  of  {u')(u')  obtained  with  three  different  k-e  closure  models 
used  in  conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  Cp(z) 
for  drag  units  ^5  and  ^6.  Model  1  (solid  line);  Model  2  (dashed  line);  Model  3  (dashed- 
dotted  line):  high-resolution  CFD  (open  circle). 
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Figure  7.  Profiles  of  ( u')(w ')  obtained  with  three  different  k-e  closure  models  used  in 
conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  Cd(z )  for 
drag  units  #1  and  #2.  Model  1  (solid  line);  Model  2  (dashed  line);  Model  3  (dashed-dotted 
line);  high- resolution  CFD  (open  circle). 
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Figure  7.  (Continued)  Profiles  of  ( u')(w ')  obtained  with  three  different  k-e  closure  models 
used  in  conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  Co{z) 
for  drag  units  $3  and  $4.  Model  1  (solid  line);  Model  2  (dashed  line);  Model  3  (dashed- 
dotted  line);  high- resolution  CFD  (open  circle). 
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Figure  7.  (Continued)  Profiles  of  ( u')(w ')  obtained  with  three  different  k-e  closure  models 
used  in  conjunction  with  the  diagnosed  values  for  the  sectional  (local)  drag  coefficient  Cd(z) 
for  drag  units  #5  and  #6.  Model  1  (solid  line):  Model  2  (dashed  line);  Model  3  (dashed- 
dotted  line);  high-resolution  CFD  (open  circle). 
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Figure  8.  Profiles  of  (u)  obtained  with  the  proposed  distributed  drag  force  model  (Model  3) 
used  in  conjunction  with  four  different  methods  for  specifying  the  drag  coefficient.  Cd{z) 
(solid  line);  CD,buik,i  (dashed  line);  CD,buik,ave  (dashed-dotted  line);  CD,buik,i  (dashed- 
dotted-dotted  line);  high-resolution  CFD  (open  circle).  Results  are  shown  for  drag  units 
#  1  and  #2. 
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Figure  8.  (Continued)  Profiles  of  (u)  obtained  with  the  proposed  distributed  drag  force 
model  (Model  3)  used  in  conjunction  with  four  different  methods  for  specifying  the  drag 
coefficient.  Cd{z)  (solid  line);  C^buik,*  (dashed  line);  C'c.bulk.ave  (dashed-dotted  line); 
C D,buik,i  (dashed-dotted-dotted  line);  high-resolution  CFD  (open  circle).  Results  are  shown 
for  drag  units  #3  and  #4. 
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Figure  8.  (Continued)  Profiles  of  (u)  obtained  with  the  proposed  distributed  drag  force 
model  (Model  3)  used  in  conjunction  with  four  different  methods  for  specifying  the  drag 
coefficient.  CD{z)  (solid  line);  C^buiM  (dashed  line);  CD.buik.ave  (dashed-dotted  line); 
CD.buik.i  (dashed-dotted-dotted  line);  high- resolution  CFD  (open  circle).  Results  are  shown 
for  drag  units  #5  and  #6. 


40 


DRDC  Suffield  TR  2004-169 


Figure  9.  Profiles  of  ( u')(u ')  obtained  with  the  proposed  distributed  drag  force  model 
(Model  3)  used  in  conjunction  with  four  different  methods  for  specifying  the  drag  coeffi¬ 
cient.  Cd(z)  (solid  line);  Co.buik.i  (dashed  line);  Cd, bulk, ave  (dashed-dotted  line);  Co.bulk.i 
(dashed-dotted-dotted  line);  high-resolution  CFD  (open  circle).  Results  are  shown  for  drag 
units  #1  and  #2. 
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Figure  9.  (Continued)  Profiles  of  ( it')  (it ')  obtained  with  the  proposed  distributed  drag 
force  model  (Model  3)  used  in  conjunction  with  four  different  methods  for  specifying  the 
drag  coefficient.  Cd{z )  (solid  line);  Co.bulk.i  (dashed  line);  C/ybuik.ave  (dashed-dotted  line), 
Co, bulk, 1  (dashed-dotted-dotted  line);  high-resolution  CFD  (open  circle).  Results  are  shown 
for  drag  units  #3  and  #4. 
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Figure  9.  (Continued)  Profiles  of  ( u')(u ')  obtained  with  the  proposed  distributed  drag 
force  model  (Model  3)  used  in  conjunction  with  four  different  methods  for  specifying  the 
drag  coefficient.  Cp(z)  (solid  line);  Co. bulk, i  (dashed  line);  Co^uik.ave  (dashed-dotted  line); 
Co.buik.i  (dashed-dotted-dotted  line);  high- resolution  CFD  (open  circle).  Results  are  shown 
for  drag  units  #5  and  #6. 
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Figure  10.  Profiles  of  ( u')(w ')  obtained  with  the  proposed  distributed  drag  force  model 
(Model  3)  used  in  conjunction  with  four  different  methods  for  specifying  the  drag  coeffi¬ 
cient.  Cd{z)  (solid  line);  Cfytmik.i  (dashed  line);  Co.buik.ave  (dashed-dotted  line),  Cx>,buik,i 
(dashed-dotted-dotted  line);  high-resolution  CFD  (open  circle).  Results  are  shown  for  drag 
units  #1  and  #2. 
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Figure  10.  (Continued)  Profiles  of  ( u')(w ')  obtained  with  the  proposed  distributed  drag 
force  model  (Model  3)  used  in  conjunction  with  four  different  methods  for  specifying  the 
drag  coefficient.  Cd(z )  (solid  line);  C^buiM  (dashed  line);  Cfybulk.ave  (dashed-dotted  line); 
Cd, bulk.  1  (dashed-dotted-dotted  line);  high- resolution  CFD  (open  circle).  Results  are  shown 
for  drag  units  #3  and  #4. 
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Figure  10.  (Continued)  Profiles  of  ( u')(w ')  obtained  with  the  proposed  distributed  drag 
force  model  (Model  3)  used  in  conjunction  with  four  different  methods  for  specifying  the 
drag  coefficient.  Cp(z)  (solid  line);  Cx>,buik,i  (dashed  line);  Cjytmik.ave  (dashed-dotted  line), 
Cx>,bulk,i  (dashed-dotted-dotted  line);  high-resolution  CFD  (open  circle).  Results  are  shown 
for  drag  units  #5  and  # 6. 
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array  of  sharp-edged  cubes  with  a  plan  area  density  of  0.25.  However,  it  should  be  emphasized  that  the 
high-resolution  CFD  flow  simulations  have  been  validated  with  wind  tunnel  experiments,  and  after  these 
validations  the  model  can  be  used  to  diagnose  the  spatially-averaged  velocity  statistics  required  for  the 
validation  of  the  distributed  drag  force  model,  it  was  found  that  the  model  predictions  for  mean  wind  speed 
and  turbulence  in  the  building  array  were  not  sensitive  to  the  differing  treatments  of  the  source  and  sink 
terms  in  the  at  -  and  E-equations,  implying  that  the  high-order  approximations  of  these  source/sink  terms 
did  not  offer  any  predictive  advantage.  A  possible  explanation  for  this  is  the  utilization  of  the  Boussinesq 
linear  stress-strain  constitutive  relation  within  the  k-z  modelling  framework,  whose  implicit  omission  of  any 
anisotropic  eddy-viscosity  effects  renders  it  incapable  of  predicting  any  strong  anisotropy  of  the  turbulence 
field  that  might  exist  in  the  building  array.  Four  different  methods  for  diagnosis  of  the  drag  coefficient  CD 
for  the  aligned  cube  array,  required  for  the  volumetric  drag  force  representation  of  the  cubes,  are 
investigated  here. 
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